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ABSTRACT 


We  investigated  an  inverse  imaging  problem  by  applying  the  Rytov 
approximation  for  Through-the-Wall  Imaging  using  THz  waves.  In  the  beginning,  we 
studied  some  properties  of  THz  waves  and  the  physical  conditions  for  THz  imaging  in 
matter.  Then  we  started  with  Maxwell’s  equations  to  derive  a  model  for  the  transmission 
of  Green’s  functions  and  used  a  Lippman- Schwinger  integral  equation  and  Rytov 
approximation  to  predict  the  scattered  field.  We  applied  the  L-curve  method  for  the 
selection  of  regularization  parameters,  and  then  presented  the  reconstruction  algorithm 
and  illustrated  the  result  with  numerical  simulations.  We  also  compared  this  result  to  the 
one  obtained  by  the  Bom  approximation. 
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I.  INTRODUCTION 


The  amazing  ability  of  terahertz  radiation  to  penetrate  through  materials  allows 
for  the  nondestructive,  noninvasive  inspection  of  mail  and  packages  in  post  offices,  and 
luggage  at  airport  security  check  points  [1]. 

The  term  terahertz  (THz)  refers  to  radiation  whose  spectrum  lies  between  the 
infrared  and  microwave  bands  with  frequencies  ranging  from  0.1  to  10  THz,  where 
1THz=1012Hz  or,  in  units  of  wavelength,  A  =3  mm  to  30  jum .  Until  recently,  this  “THz 
gap”  has  been  almost  inaccessible  due  to  the  lack  of  efficient  sources  and  detectors  in  this 
region  [2].  But  recently,  sources  and  detectors  for  THz  frequencies  have  become 
available  and  easier  to  use.  Therefore,  many  applications  of  terahertz  radiation,  such  as 
imaging  and  communications,  are  now  appearing.  Among  these  applications,  a  very 
important  application  (and  the  focus  of  much  research)  is  imaging  using  THz  radiation 
[3].  The  study  of  this  imaging  is  an  interesting  and  attractive  area  for  researchers. 

A.  BACKGROUND 

Various  THz  imaging  techniques  have  been  developed  since  Hu  and  Nuss  first 
demonstrated  THz  imaging  in  1995  [4].  THz  imaging  applications  are  the  focus  of 
constantly  growing  interest.  Because  THz  waves  can  provide  great  spatial  resolution  and 
can  penetrate  most  dielectric  materials,  such  as  plastics,  paper,  and  wood,  there  are  many 
applications  of  THz  imaging.  A  classical  application  of  THz  wave  based  remote  sensing 
is  content  inspection  of  packages  and  envelopes.  Recently,  THz  imaging  has  been  used 
for  medical  and  biological  applications  because  it  provides  high  quality  images.  Another 
application  is  real-time  nondestructive  testing,  which  includes  testing  for  defects  in 
plastic  food  and  medicine  packages.  Also,  there  are  a  number  of  security  applications  for 
THz  imaging:  for  example,  Through- the- Wall  Imaging  (TWI),  luggage  inspection,  and 
the  detection,  from  a  distance,  of  weapons  or  explosives  hidden  under  clothing  or 
briefcases. 
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B.  MOTIVATION 

In  the  decade  between  1996  and  2005,  more  than  a  half  million  servicemen  were 
injured  or  killed  in  situations  involving  barricaded  offenders,  hostage-taking,  and  high- 
risk  entry  [5],  Through-the-Wall  Imaging  (TWI)  has  been  studied  for  saving  soldier’s  and 
officer’s  lives.  TWI  technology  can  help  soldiers  to  detect  enemies  hiding  in  a  building, 
and,  for  this  reason,  TWI  technology  needs  to  evolve  and  improve. 

The  armed  forces  personnel  or  police  officers  who  experience  operations  on 
urbanized  terrain  may  want  equipment  that  can  detect,  identify,  and  localize  enemies 
through  a  wall,  and  operate  day  and  night  in  all  weather  conditions. 

Most  imaging  systems  use  bi-static  or  multi-static  transmission,  where  the 
transmitter  is  on  one  end,  and  the  receiver  is  on  its  opposite  end.  Mono-static  systems,  on 
the  other  hand,  use  transmitters  and  receivers  that  are  in  the  same  location.  Mono-static 
imaging  systems  rely  on  the  back-scattered  field.  Through- the- wall  imaging  is  more 
difficult  and  challenging  than  any  other  back  scattering  method  because  of  the  general 
variability  of  wall  construction.  In  this  thesis,  we  will  concentrate  on  mono-static  imaging 
techniques. 

There  are  two  basic  problems  in  ray-scattering.  The  first  one  is  the  direct 
scattering  problem,  which  predicts  how  the  rays  scatter  from  a  known  object.  The  second 
is  the  inverse  scattering  problem,  which  attempts  to  recover  the  original  object  from  the 
blurred  and  corrupted  measurements  of  the  scattered  field.  The  inverse  scattering  problem 
is  one  of  determining  the  characteristics  of  an  object  from  measurement  data  of  radiation 
or  particles  scattered  from  the  object.  Also,  most  inverse  problems  are  ill-posed.  X-ray, 
Tomography,  and  ultra  sound  imaging  are  good  examples  of  ill-posed  inverse  problems. 

Let’s  look  at  the  difference  between  the  direct  problem  and  the  inverse  problem. 
The  forward  problem  can  be  expressed  as 

D  =  A(f)  (1.1) 


in  which,  we  denote  the  data  by  D  e  □  M  ,  the  object  by  f  e  N  and  the  operator  acting 
between  two  spaces  containing  D  and  /  by  A .  ( N  and  M  denote  the  dimensions  of, 
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respectively,  the  object  and  data  spaces  and,  in  general,  N  *  M  .)  So  the  direct  problem 
predicts  the  data  (D)  from  the  object  (/).  On  the  other  hand,  the  inverse  problem  is 
expressed  as 

/  =  A\D)  (1.2) 

and  solves  for  the  object  (/)  from  the  data  (D).  This  equation  looks  easy  to  solve,  but 
real-world  inverse  problems  are  typically  much  more  complicated. 

Now  let’s  look  at  the  inverse  problem  a  little  bit  deeper.  If  we  want  to  accurately 
map  objects  to  data,  we  must  also  consider  the  measurement  errors.  The  difference 
between  the  actual  measurement  data  ( D  )  and  the  real  data  ( D )  is  denoted  by  noise 
n=D-D ,  and  the  scattering  problem  can  be  more  accurately  expressed  by 

D  =  A(f)  +  n  (1.3) 

The  object  and  data  functions  are  conveniently  considered  as  belonging  to  Hilbert 
spaces.  In  equation  (1.3),  A  is  an  operator  from  one  Hilbert  space  to  the  other.  Hilbert 
space  allows  simple  geometric  concepts  like  projection  onto  fewer  dimensional  spaces 
and  loss  of  information. 

Inverse  problems  are  almost  always  ill-posed  problems.  J.  Hadamard  fonnalized 
the  concept  of  well-posedness  for  equations  of  the  fonn  (1.3)  [6].  Any  equation  that  does 
not  satisfy  all  three  conditions  is  called  ill-posed,  and  a  method  for  solving  this  problem 
approximately  is  called  a  regularization  method. 


1 

Existence 

There  exists  a  solution  to  the  problem 

2 

Uniqueness 

There  is  at  most  one  solution  to  the  problem 

3 

Stability 

The  solution  depends  continuously  on  the  data 

Table  1.  Conditions  for  Well-posed  Problems. 
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c. 


PROPERTIES  OF  TERAHERTZ  FREQUENCIES 


The  term  “Terahertz”  is  applied  to  broadband  pulsed  electromagnetic  radiation 
with  a  spectrum  that  falls  between  the  infrared  and  microwave  bands  ranging  from  0.1  to 
10  THz  [7].  The  wavelength  in  this  region  ranges  from  3  mm  to  30  jum  (as  shown  in 
figure  1).  These  are  short  wavelengths,  and  they  have  high  resolution  compared  to  other 
radio  frequencies  [2],  There  are  several  distinct  advantages  of  THz  waves.  The  first  is 
that  most  chemical  substances  possess  characteristic  absorption  features  in  the  THz  range. 
This  property  enables  us  to  detect  illegal  chemicals,  even  when  they  are  sealed  inside  a 
packet  or  clothing.  The  second  advantage  is  material  transmission  properties.  Most  non¬ 
conducting  materials  such  as  plastics,  paper,  wood,  and  ceramics  are  totally  or  partially 
transparent  in  the  THz  range.  The  third  advantage  is  that  the  radiation  is  non-ionizing  and 
has  relatively  low  energies  compared  to  X-rays;  therefore,  we  don’t  have  to  be  concerned 
about  the  safety  of  THz  imaging. 

In  this  thesis,  we  will  use  THz  pulses  for  imaging.  Figure  2(a)  shows  examples  of 
THz  pulses  transmitted  through  nylon.  When  a  THz  pulse  interacts  with  an  object,  it 
experiences  a  delay,  an  attenuation,  and  a  broadening,  as  the  different  component 
frequencies  are  phase  shifted,  absorbed,  reflected,  and  scattered.  Figure  2(b)  shows  the 
power  spectrum  that  the  pulse  frequency  content  is  centered  around  ITHz  (note  that  the 
higher  frequency  content  is  attenuated  more  than  the  lower  frequencies  by  nylon). 


Figure  1.  Electromagnetic  spectrum[From  7] 
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Figure  2.  Examples  of  terahertz  pulses  in  (a)  Time  Domain  and  (b)  Magnitude  of 

Fourier  Coefficients  [From  4] 

For  the  convenience  of  modeling,  we  assume  that  the  wall  that  we  want  to  see 
through  is  a  linear,  homogeneous,  and  isotropic  dielectric.  We  know  THz  frequencies 
have  good  penetration  for  dry  materials,  but  a  spectral  cut-off  above  3  THz  will  be  caused 
by  absorption  from  water  vapor  in  the  THz  beam  path  [8].  So  we  are  not  generally 
interested  in  frequencies  over  3  THz.  We  also  need  to  consider  the  object  that  we  want  to 
image.  The  reflection  from  the  object  depends  on  the  object  properties.  For  example, 
metallic  materials  reflect  perfectly,  but  skin  or  water  absorbs  THz  rays,  so  we  can’t  get 
significant  reflected  pulses.  Therefore,  we  should  be  careful  to  choose  frequencies  that 
balance  between  the  penetration  of  medium  and  the  reflection  of  the  object. 

D.  THESIS  ORGANIZATION 

This  thesis  uses  Born  and  Rytov  approximations  to  model  measured  scattered 
fields  from  objects  located  behind  walls. 
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Chapter  II  introduces  the  concept  of  inverse  problems,  including  ill-posed 
problems,  singular  value  decomposition,  and  regularization  methods.  The  Born  and 
Rytov  approximations  are  discussed  as  well. 

Chapter  III  provides  some  properties  of  electromagnetic  waves  in  media  and 
restricts  our  model  system  to  a  real  situation.  A  scalar  theory  is  introduced  for  use  in  the 
succeeding  chapters. 

Chapter  IV  develops  the  transfer  function  of  the  media.  This  chapter  describes  the 
mathematical  approach  and  derives  Green’s  functions  in  different  regions. 

Chapter  V  examines  reflectives  of  media  and  evanescent  waves.  It  also  examines 
materials  with  reflectivity. 

Chapter  VI  derives  a  wave  equation  in  the  specified  medium  and  the  scattered 
fields  by  applying  the  Born  and  Rytov  approximations. 

Chapter  VII  provides  algorithms  and  simulations  (with  noise)  by  applying 
Tikhonov  regularization  and  Truncated  Singular  Value  Decomposition  regularization 
schemes. 

Finally,  Chapter  VIII  gives  a  final  summary  conclusion,  comparing  the  simulation 
results  and  providing  suggestions  for  future  work. 
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II.  THE  INVERSE  PROBLEM 


This  chapter  will  introduce  ill-posed  problems  and  regularization  solution 
methods  more  deeply.  Regularization  methods  are  used  to  mitigate  noise  effects  in  the 
data  or  the  operator.  Three  widely  used  regularization  techniques  are  Tikhonov 
regularization.  Truncated  Singular  Value  Decomposition,  and  regularizing  iterative 
methods.  In  this  thesis,  we  will  use  the  Tikhonov  regularization  method  and  the 
Truncated  SVD  method.  We  will  confine  ourselves  to  linear  equations  with  compact 
operators  in  Hilbert  space  and  use  the  Singular  Value  Decomposition  for  our 
reconstruction  algorithm  design. 

A.  REVISIT  ILL-POSED  PROBLEMS 

We  want  to  get  images  using  electromagnetic  frequencies.  So  there  is  no  way  to 
avoid  the  inverse  problem.  The  very  first  equation  we  need  is  a  Lippman-Schwinger 
equation  [9]. 

Vscatti*)  =  JJLG^(x’xV(V)(^„c(V)  +  ii/scatt(x'))d:'x'  (2. 1) 

We  denote: 

•  Gk{x,x  )  is  Green’s  function  which  satisfies  the  wave  equation  with  an 
impulsive  source  term 

•  p{x  )  is  a  source  factor  detennined  by  the  index  of  refraction  of  the 
scatterer 

•  y/{x)  is  an  electric  field,  y/inc(x)  is  the  incident  field,  and  ipscall(x)  is 
the  scattered  field 

Equation  (2.1)  is  derived  from  the  reduced  scalar  wave  equation 

V2y/(x)  +  k2n2(x)y/(x)  =  0  (2.2) 

where  n(x)  is  the  index  of  refraction  and  k  is  the  wave  number. 
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In  this  thesis,  we  ignore  multiple  scattering.  Because  we  assume  multiple 
scattering  terms  are  very  small  in  comparison  to  primary  scattering,  terms  we  can  rewrite 
equation  (2.1)  as 

Yscatt  O)  =  \\\DGk{x,x')p{x')¥inc(x')d'x'  (2.3) 

To  get  images  of  objects,  we  have  to  estimate  p{x)  (which  defines  the  “image”), 
but  there  are  some  problems.  The  main  problem  is  that  we  can  never  obtain  data  for  all 
points  (x)  and  at  all  frequencies  (co) .  The  other  problem  is  that  the  measurements  are 
always  accompanied  by  noise.  If  we  consider  these  problems,  we  can  make  a  more 
complete  measurement  model  and  express  equation  (2.3)  in  matrix  form.  To  perform  this 
discritization,  we  break  the  integral  in  equation  (2.3)  interval  x  e  [a  b\  into  M  parts  and 
ie[c  d]  into  N  parts.  Then  the  integral  equation  becomes 

rb 

K(x.  5)  /  ( s)ds  =  d(x.)  a<s<b,  c  <  x.  <  x2  <  •  •  •  <  xN  <  d  (2-4) 

J  a  ’ 

where  operator  K  is  defined  as  the  kernel  which  operates  on  the  object  / el  to 
produce  the  measurement  d  eY  .  We  have  simplified  equation  (2.3)  into  a  one 
dimensional  case  and  changed  notations  to  better  understand  this  measurement  system. 

Where, 

•  The  measurement  at  position  xt  :  <//s<;a„(xi)  =>  d(x.) 

•  The  “kernel”  :  Gj  (x, x')  =>  K(x;.,s) 

•  The  “object”  function  :  p{x  )  =>  f(s ) 

Then  we  can  rewrite  equation  (2.4)  in  matrix  form  (with  noise  included)as 

d  =K f  +  n  (2.5) 

where 
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d(x j) 

K(xj  a  +  As) 

K(xj  a  +  2As) 

••  K(xja-t-MAy) 

d  = 

d(x 2) 

K  = 

K(x2a  +  As) 

K(x2  a  +  2A  s)  ■ 

•  •  K(x,  a  +  MAs ) 

_d(xN)_ 

K(xw  a  +  Ay) 

K(xa,  a  +  2Ay)  • 

••  K(xw  a  +  MAs) 

f(a  +  As) 

nix  j) 

f  = 

f(a  +  2As) 

,  and  n  = 

n(x2) 

f{a  +  MAs) 

_n(xN)_ 

and  d  e  and  /  e  Mm  are  vectors. 

Let’s  take  a  close  look  at  the  dimension  of  the  measurement  space  and  the  object 
space.  The  measurement  space  has  dimension  N ,  and  the  object  space  has  dimension  M  . 
This  reflects  the  number  N  of  the  discrete  measurements  and  the  number  N  defining  the 
resolution  with  which  /  can  be  estimated.  Normally  the  object  dimension  M  is  greater 
than  the  measurement  of  dimension  N .  Now,  equation  (2.5)  can  be  used  to  get  /  by 
defining  an  inverse  K  .  However,  the  noise  can  alter  the  measurement  space  and  lead  to 
an  unstable  estimate  of  /  =  K d'd .  Equation  (2.5)  is  also  often  an  over-detennined  system, 
which  means  there  are  more  equations  than  unknowns,  and  we  can  estimate  /  by  least- 
squares.  Therefore,  we  have  to  figure  out  an  over-determined  system  to  get  solutions 
close  to  the  real  (discrete)  values.  We  will  deal  with  some  methods  to  get  stable  solutions 
in  section  D. 


B.  LEAST  SQUARES  SOLUTION  BY  SINGULAR  VALUE 
DECOMPOSITION 


We  consider  d  to  be  a  vector  with  components  di  =  d(xi).  Normally,  the  vector 
of  measurements  d  in  equation  (2.5)  will  not  lie  exactly  in  the  measurement  space 
Y  (because  of  noise).  Because  d  <eY  only  if  the  noise  vector/?  e  Y  and  K f  e  Y .  But  n  is 
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Figure  3.  Hilbert  space  and  measurement  space 

independent  of  the  operator  K  ,  and  the  object  vector  /  may  have  components  in 
measurement  space  or  in  null  space.  So  a  possible  “best  object”  f  will  minimize  the 
distance  \\d  — K/'|| .  We  denote  the  least  squares  solution  to  equation  as 

/  =  arg  min \\d  - K/||  (2.7) 

From  Figure  3  we  can  get  the  solution  to  equation  (2.7).  The  vector  K f  lies  in 
measurement  space  Y  .  When  K f  =  P  ,  K f  gets  the  closest  distance  to  d  .  The  vector 
d  -P  =  d  -Kf  is  orthogonal  to  the  measurement  space  Y  and  all  vectors  of  the  fonn  Kv 
in  measurement  space  Y .  Then  we  can  get  the  normal  equations  by  [10] 

KTKf  =  KTd  (2.8) 

/  is  the  least  square  solution  to  equation  (2.8).  Since  the  matrix  KrK  e  □  M/M  is  square 
we  can  write 

f  =  (KTK)-1KTd  (2.9) 

/  is  exactly  what  we  want.  In  order  to  get  / ,  we  have  to  compute  (KrK)~'Kr  .  But  this 
can  be  painful  to  invert  with  large  N  and  M  .  So  we  want  to  use  another  method, 
singular  value  decomposition,  to  get  / . 
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We  would  like  to  employ  simple  linear  algebra  to  build  matrices  of  the  form  (KrK)  1 K.7  . 
Consider  a  matrix  A  that  is  square  (M  x  M )  with  e}  an  eigenvector,  2,  an  eigenvalue. 
We  can  transfonn  the  matrix^  in  tenns  of  the  matrices  e;ej ,  formed  from  the  set  of 
vectors  {e,.} ,  using  the  Spectral  Theorem  which  states  that  A  can  be  written  with  as 

M 

A  =  UAU ~l  =  UAUr  =  YAietef  (2.  i°) 

;=i 

where  U  is  the  matrix  with  the  eigenvectors  of  A  on  its  columns  and  A  is  the  matrix 
with  the  eigenvalues  of  A  along  its  diagonal. 

We  assume  {e,}  is  an  orthonomal  basis  of  □ M  ,  and  so  we  can  represent  an 
M  dimensional  vector  v  in  terms  of  e: 

M  M 

y  =  Y  (efy>i  =  £  yfit  (2- 1 1) 

i= 1  i=l 

where  y.  =  ef  y  is  the  component  of  the  vector  y  in  the  basis.  Consequently, 

M  MM 

x  =  Ay  =  Y  ACCT  =  Y  AAC  =  Y  AC  (2- 12) 

i= 1  i=l  i= 1 


This  is  representation  of  x  =  Ay  in  terms  of  the  basis  ( e, } . 


Let’s  adjust  this  result  for  the  matrix  (KrK)  1 K7  .  Suppose  are  the  eigenvectors  of 
KrK .  Then 


(K^K)-1  =  (UAU^y1  =  UA~lUr 


(2.13) 


Multiply  by  Kr  on  the  right  of  both  sides 


11 


(KrKr*Kr  =  UA~1UTKT  =UA-\KU)t 


f 

X 

UA-1 

V 

1 

...  ^ 

A1 

...  ^ 

••  K*  1 

J 

1 

q 

> 

0 

1 - 

o 

1 

■  •  o 

CTj/zf 

•  •  o 

...  eT 

1 

O  • 

•  o 

r 

b" 

...  eT 

L  M 

D  VT 


=  udvt 


We  define  cr.  =  Ke,.  and  cr.  =  ex 


jo-r'Ke,.  if  cri  *  0 
[0  if  cr  =  0 


Now  multiply  both  sides  of  equation  (2.9)  by  K  so  that 

K7  =  K(KrK)  1KV=P 


(2.14) 


(2.15) 


P  is  the  projection  of  the  measurement  vector  onto  measurement  space.  Therefore, 
K(KrK)  1 K7  is  a  projection  matrix  mapping  d  to  its  components  in  measurement  space. 
We  will  use  this  projection  matrix  to  get  more  accurate  data  using  regularization 
techniques  in  the  following  section. 


C.  TRUNCATED  SINGULAR  VALUE  DECOMPOSITION 


To  obtain  a  better  estimate  of  the  least  square  solution,  the  truncated  singular 
value  decomposition  solution  is  often  used.  This  method  truncates  the  SVD 
representation  to  reduce  the  effect  of  noise  contamination.  Just  select  one  small  threshold 
value  JC  and  compare  this  with  singular  values  f .  If  /C  >  f . ,  then  we  replace  the  1 / 

in  D  by  0.  In  this  thesis,  using  the  plot  of  singular  values,  we  will  find  the  critical  point 
where  the  singular  values  change  steeply  and  will  truncate  at  that  point  [11]  (see  Figure 
21). 


12 


D.  REGULARIZATION  METHODS 


Since  our  model  is  linear  and  ill-posed,  the  inverse  image  is  under-determined, 
and  we  will  get  small  singular  values  of  K.  The  problem  is  that  K  depends  on  our  model 
of  the  measurement  process  and  is  not  completely  known.  So  we  get  singular  values  with 
a  slight  imprecision.  Now  we  introduce  regularization  methods  to  overcome  the  problems 
associated  with  the  small  singular  values  that  lead  to  reconstruction  instability. 
Regularization  is  important  for  solving  inverse  problems,  because  the  output  of  least 
squares  is  affected  by  data  and  rounding  errors.  Therefore  we  would  like  to  introduce 
regularization  methods  to  reduce  these  errors.  But  there  is  a  trade-off  between  the 
regularized  output  and  the  original  sets  of  data.  If  the  regularization  is  too  crude,  the 
regularized  solution  does  not  fit  the  given  data,  and  the  residual  error  ||r||  =  \d  -K/"||  is 

unnecessarily  large.  If  the  regularization  is  too  small,  the  fit  will  be  good,  but  data  noise 
effects  will  be  larger. 

There  are  also  some  regularization  methods  that  employ  operator  or  data 
correction  [12], 

1,  Regularization  by  Operator  Correction 

Let’s  model  the  problem  by  mapping  A  as 

Ax  =  d,  xeX,  d  eY  (2-16) 

The  equation  (2.16)  is  said  to  be  an  ill-posed  problem,  if  it  does  not  meet  one  of  the 
conditions  in  table  I.  If  we  want  to  transform  an  ill-posed  problem  into  a  well-posed 
problem  by  approximation  of  the  equation,  then  we  have  to  choose  a  mapping  method. 

The  regularization  of  equation  (2.16)  consists  of  a  substitution  of  the  operator  A 
by  A  .  Let  A-  X  — »  7  be  an  operator  such  that 

Ax  =  d  (2.17) 

This  is  the  idea  behind  the  well  known  Tikhonov  regularization. 
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2.  Regularization  by  Data  Correction 


Let  K-Y  — »  A  be  a  continuous  operator  such  that 

then 

Ax  =  Kd,  xeX,  deY  (2.18) 

is  well-posed  and  called  as  regularization  by  data  correction. 

3.  Regularization  by  Data  and  Operator  Correction 

The  last  method  of  regularizations  is  a  combination  of  both  previous  methods.  Let 
A-  X  — »  7  be  an  operator  and  K  -  Y  —>  Y  a  continuous  mapping  such  that 

Ax  =  Kd,  xeX,  deY  (2.19) 

Tikhonov  regularization  is  also  a  well  known  method  of  this  kind  of  regularization. 

E.  L-CURVE  METHOD 

Since  the  quality  of  the  regularized  solution  is  controlled  by  the  regularization 
parameter  for  example,  the  threshold  parameter  /C  used  in  truncated  SVD  regularization 
[13],  we  have  to  choose  the  optimal  parameter  to  get  a  “best”  image.  There  are  various 
methods  for  the  selection  of  regularization  parameters,  including  the  Discrepancy 
Principle,  Generalized  Cross  Validation,  and  L-curve.  The  L-curve  method  applies  a  log- 
log  plot  of  the  regularized  solution  against  the  squared  norm  of  the  regularized  residual 
for  a  range  of  values  of  regularization  parameters  [13]  and  selects  the  corner  as  the  point 
of  maximum  curvature  in  the  L-curve  plot  (see  Figure  4). 

There  are  two  main  approaches  used  in  the  L-curve  method.  The  first  approach 
considers  both  the  residual  nonn  and  the  solution  norm.  The  second  approach  considers 
the  value  of  the  maximum  curvature.  Like  other  methods,  the  L-curve  method  has  its 
merits  and  limitations.  The  merits  of  L-curve  are  that  it’s  robust  and  can  treat 
perturbation  caused  by  correlation  noise.  On  the  other  hand,  the  L-curve  method  is 


14 


limited,  because  it  is  not  convergent  [14],  i.e.,  since  singular  values  are  getting  close  to  0 
but  not  exactly  0  as  increasing  index  of  singular  values,  we  can’t  get  the  value  which 
satisfies  the  norm  is  0. 


Figure  4.  L-curve  for  Tikhonov  Regularization  [From  14] 


F.  APPROXIMATIONS  TO  THE  WAVE  EQUATION 

In  the  case  of  electromagnetic  scattering  problems,  we  need  to  approximate  the 
unknown  internal  field  by  a  known  distribution  in  order  to  linearize  the  problem.  The 
Born  and  Rytov  approximations  are  widely  used  to  simplify  electromagnetic  scattering 
problems.  In  this  thesis  we  will  derive  reconstruction  formulas  under  the  Rytov 
approximation  and  compare  with  the  Bom  approximation. 


1.  Born  Approximation 


The  first  Born  approximation  was  introduced  by  Max  Born  in  1925,  in  order  to 
solve  problems  concerning  the  scattering  of  atomic  particles  [7].  The  Born  approximation 
is  obtained  by  expressing  the  total  wave  field  as  the  sum  of  incident  field  plus  scattered 
field: 


y(r)  =  y/i(r)  +  y/s(r) 


(2.20) 
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This  approximation  assumes  that  scattered  fields  inside  scatterers  are  negligible 
compared  to  the  background  and  the  normal  fields.  The  first-order  Bom  approximation  is 
good  only  if  the  scattered  field  is  much  smaller  than  the  incident  field. 

2.  Rytov  Approximation 

The  Rytov  approximation  was  discovered  in  1937  by  Rytov,  who  was  analyzing 
the  diffraction  of  light  by  sound  waves.  This  approximation  was  later  applied  by 
Obukhov  to  describe  the  propagation  of  electromagnetic  waves  in  random  media  [15]. 
The  Rytov  approximation  is  based  on  expressing  the  total  field  as  a  complex  phase, 
which  is  the  sum  of  incident  phases  plus  scattered  phases: 

^(r)  =  /r)=e*(rW(r)  (2.21) 

The  condition  for  the  applicability  of  this  approximation  is  that  the  phase  of  the  scattered 
field  varies  slowly,  relative  to  one  wavelength.  So  the  size  of  target  is  less  critical  in  the 
Rytov  approximation. 


16 


III.  ELECTROMAGNETIC  WAVES  IN  MEDIA 


It  is  difficult  to  explain  electric  fields  and  magnetic  fields,  which  are  cross- 
coupled,  without  the  Maxwell’s  equations.  Maxwell’s  equations  were  established  by 
James  Clerk  Maxwell  by  1864  and  experimentally  verified  by  Heinrich  Hertz  in  1888 
[16].  Since  then  these  equations  have  been  applied  in  a  variety  of  areas  such  as  optics, 
microwaves,  antennas,  radar,  and  communications.  In  this  chapter  we  will  investigate 
Maxwell’s  equations  in  media  and  derive  the  associated  wave  equations. 


A.  MAXWELL’S  EQUATIONS 


Maxwell’s  equations  are  comprised  of  Gauss’  law  for  the  electric  field,  the 
observation  that  there  are  no  magnetic  monopoles,  Ampere’s  law,  including  the 
dD 

displacement  current  — ,  and  Faraday’s  law  of  induction  in  differential  form.  They  can 

dt 

be  written  as 


W-D  =  p 
V-B  =  0 


VxE 


dB 

dt 


V  xH  =  J  + 


dD 

dt 


(3.1) 


where  D  and  B  are  the  electric  and  magnetic  flux  intensities,  and  E  and//  are  the 
electric  and  magnetic  field  intensities,  J  is  electric  current  density,  and  p  is  volume 
charge  density.  In  reality,  there  are  two  more  things  we  have  to  consider  for  Maxwell’s 
equations.  These  are  P ,  and  M  ,  called  induced  polarization  and  magnetization 
respectively.  The  associated  constitutive  relations  are  expressed  as 


D  —  +  P 

(3.2) 

B  =  jU0(H  +  M) 

(3.3) 
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Using  equations  (3.2)  and  (3.3),  we  can  rewrite  Maxwell’s  equations  in  terms  of  E  and B  . 


V-E  =  —(p-V-P) 

eo 

V-B  =  0 


VxE 


dB 

dt 


„  „  dE 
VxB  =  p0e0  — +  ju0 
dt 


dP 

J  +  —  +  VxM 
dt 


(3.4) 


where  ju0  and  e0  are  the  permeability  and  pennittivity  in  free  space,  respectively.  The 
speed  of  light  (c)  and  characteristic  impedance  (^0)  of  a  vacuum  are  expressed  using  ju0 
and  e0  as 


c  = 


(3.5) 


1.  Simplified  Model 


Electromagnetic  fields  behave  differently  in  the  presence  of  media.  For  Through 
the  Wall  Imaging  using  electromagnetic  waves,  we  will  confine  our  attention  to  materials 
obeying  a  linear  model.  Materials  used  for  constructing  walls  include  wood,  concrete, 
and  iron-beams.  These  materials  typically  have  anisotropic  permittivity.  Anisotropic  is  an 
inherent  property  of  the  atomic  and  molecular  structure  of  the  dielectric  [17].  So  in 
anisotropic  materials,  permittivity  depends  on  its  direction  and  the  linear  constitutive 
relation  can  be  written  in  matrix  form  as 


X  1 

e 

XX 

ev; 

XI 

Dy\ 

= 

e„ 

% 

evz 

Ey\ 

(3.6) 

P:\ 

A 

% 

e,. 

X  J 

In  general,  the  field  vector  E  is  no  longer  parallel  to  D  .  But  we  can  simplify  our  model 
by  using  the  principal  axis  system:  We  consider  homogeneous  and  linear  materials  so  that 
the  permittivity  is  independent  of  position.  Then  the  permittivity  tensor  can  be  expressed 
with  six  elements  as 
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tensor (q )  = 


xx  xy  xz 


xy  yy  yz 


xz  yz  zz 


(3.7) 


Since  this  matrix  is  symmetric,  its  eigenvectors  are  orthogonal.  If  we  choose  a 
coordinate  system  which  is  aligned  to  these  eigenvectors,  we  can  rewrite  the  permittivity 
tensor  as  a  diagonal  tensor  determined  by  its  eigenvalues.  Then  equation  (3.7)  becomes 


tensor  (e)  = 


e  0 

XX 

0  e, 
0 


yy 

0 


0 

0 

e 


(3.8) 


Let’s  investigate  solutions  to  the  Maxwell’s  equations  for  source  free 
homogeneous  isotropic  media.  Here,  J  =  p  =  0  and  Maxwell’s  equations  (3.4)  become 


V-E  =  0,  VxE  +  —  =  0, 
dt 

r)E 

V-B  =  0,  VxB-ue —  =  0 

8t 


(3.9) 


To  derive  the  wave  equation  we  first  take  the  curl  of  the  second  of  equation  (3.9) 

Vx(Vxf)  +  Vx  —  =  0 
dt 


(3.10) 


Vx(Vx£)  +  — (VxB)  =  0 
dt 


(3.11) 


and  use  the  identity  of  V  x  (VxE)  =  V(V  -E)~  V2E  =  -V2E,  since  we  assume  a  charge 
free  and  homogeneous  region  that  means  V  •  E  =  0 ,  and  we  obtain 


t  d 

-V‘T  +  -(Vx5)  =  0 
dt 


(3.12) 


Now  take  the  partial  derivative  with  respect  to  t  of  the  forth  of  equation  (3.9) 


vy  dB  dlE  A 

V  x - ue — r-  =  0 

dt  dt2 


i,(VxB)-fiew 


(3.13) 


and  substitute  this  into  equation  (3.12)  to  obtain 
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(3.14) 


V2E  -  jue(-^-  =  0 

dr 

The  plane  wave  solutions  to  equation  (3.14)  can  be  written  as 

E(r,t)  =  E0el(lp-m)  (3.15) 


Substituting  equation  (3.15)  into  equation  (3.14)  we  find  the  following  equation 

k2  =  orjue  (3.16) 


This  is  called  as  the  dispersion  relation.  We  can  easily  find  the  phase  velocity  as 


a) /k  =  1 Jyfjue 


(3.17) 


For  conductive  materials  with  permittivity  e+  i  a / co ,  the  wave  number  is  complex  and 
written  as  [16] 


k  =  kR+  ik, 


.  a 
/  — 
coe 


(3.18) 


Then  the  penetration  depth  is 


2 

COjUCT 


(3.19) 


For  a  wave  propagating  in  the  +z  direction,  we  have  z  dependence  given  by 

g-'fe  _  e-'kRz~'kiz 


(3.20) 


This  wave  attenuates  exponentially  in  the  same  direction  that  it  propagates.  Figure 
5  shows  a  wave  propagating  in  a  lossy  medium. 

,  cr 
1<  — 

For  a  highly  conducting  medium  with  coe ,  the  penetration  depth  is  a  very 
small.  So  the  waves  will  not  be  able  to  penetrate  conductive  materials.  Therefore  we 
ignore  highly  conductive  materials  in  this  thesis. 
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Figure  5.  Wave  propagating  in  a  lossy  medium[From  18] 

2.  Properties  of  Frequency  and  Dispersion 

When  the  electric  field  is  applied  to  a  dielectric,  the  electric  field  tends  to  separate 
the  electron  from  the  positively  charged  nucleus,  and  this  creates  an  electric  dipole 
moment.  Moreover,  the  dipole  moment  contributes  to  the  electric  flux  density  so  that 

D  =  e0E  +  P  =  e0(\  +  x(co))E  =  e(co)E  (3.21) 

where  x(co)  is  susceptibility.  From  equation  (3.21)  we  can  see  the  permittivity  is 
changed  because  of  the  dipole  moment.  The  polarization  of  a  medium  with  refractive 
index  n  is 

P  =  X%E 

and  the  susceptibility  is  expressed  as  [16] 

2 

nq 

X  = - - 2  rT 

me0(co0  -co  -ly) 


(3.22) 


(3.23) 
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where  y  =  J/ ,  which  is  a  measure  of  the  rate  of  collisions  per  unit  time.  Rationalizing 
Equation  (3.23)  as 


(3.24) 


and  substituting  equation  (3.24)  into  equation  (3.21),  we  get  the  pennittivity  as 


2  t  2  2  \ 

nq  [co0  -co  I 


inq  coy 


m 


(2  2\ 

2  i  2  2 

-  + — 

(2  2\ 

2  2  2 

[C00~CD  ) 

+  co  y 

m 

[C00~CD  j 

+  co  y 

(3.25) 


When  co  is  not  close  to  co0 ,  the  imaginary  tenn  will  be  much  smaller  than  the  real  term. 

This  means  e  is  approximately  real  valued  when  co  is  not  very  near  resonance.  When  not 
close  to  a  resonance,  the  deviation  of  the  refractive  index  from  unity  is 


n-m 


(  ■  \ 

1+  m. 

2  2 

{  G)0  -  CO  J 


(3.26) 


From  equation  (3.26)  the  real  part  of  n  represents  dispersion  and  the  imaginary  part 
represents  absorption.  Figure  6  shows  us  that  around  the  resonant  frequency  co0 ,  the  real 
part  of  n  behaves  in  a  strange  manner  and  drops  rapidly  with  frequency,  and  the  material 
absorption  occurs  quite  high.  So  we  have  to  avoid  choosing  a  radiation  frequency  that  is 
in  resonance  with  the  dielectric. 


To  make  our  model  system  more  numerically  tractable,  we  need  more 
assumptions.  Dispersive  media  have  values  of  e  and  °  that  depend  on  frequency.  As 
a  result,  the  wave  velocity  is  generally  frequency-dependent  [18].  So,  dispersive  media 
have  frequency-dependent  constitutive  laws;  this  means  waves  of  different  frequencies 
propagate  with  different  velocities.  For  dispersive  media,  we  need  lots  of  information 
about  different  materials.  The  analysis  is  complex  and  unnecessary.  So  we  will  not 
consider  highly  dispersive  media  in  this  thesis. 
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Figure  6. 


Real  and  imaginary  parts  of  the  refractive  index  [From  19] 


In  general,  the  velocity  of  the  envelope  of  any  modulated  sinusoid  is  the  group 
velocity  Vg  which  is  the  speed  at  which  energy  and  information  travel  [18]. 


Vg  dco 


8k_ 

ydco  j 


(3.27) 


For  these  reasons,  our  model  assumes  that  the  media  have  low  conductivity  and  are  non- 
dispersive. 


B.  THE  SCALAR  THEORY 

In  the  previous  section,  we  derived  the  vector  wave  equation  (3.14).  Also,  we 
restricted  our  model  to  a  dielectric  medium  that  is  linear,  homogeneous,  isotropic,  and 
non-dispersive.  If  the  medium  is  just  like  that,  we  can  rewrite  the  wave  equation  for  the 
electric  field  as 


V2E 


n  d2E 

~7W 


=  o 


Identically,  we  get  the  same  equation  for  the  magnetic  field  as 


V2H 


n  o  H 


=  0 


(3.28) 


(3.29) 


Since  these  are  vector  equations,  we  will  get  6  equations  with  6  components:  Ex,E  ,E_ 
and  Hx,H  ,H_ .  But  fortunately,  we  can  make  these  equations  into  a  single  scalar  wave 
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equation  because  under  our  model  restriction,  all  vector  components  of  E  and  H  are 
identical.  A  single  wave  equation  becomes 


vyo,o  + 


n  82i//(x,t) 
c7  dt2 


=  0 


(3.30) 


where  y/(x,t)  represents  any  of  the  vector  field  components  and  t//  is  dependent  on 
position  x  and  time  t .  However,  when  the  medium  is  inhomogeneous  with  a  permittivity 
e(jc)  that  depends  on  position  x ,  then  the  wave  equation  is  not  the  same  as  before.  The 
wave  equation  is  derived  by  the  following  procedures.  First  we  consider  the  first  line  of 
equation  (3.1)  and  equation  (3.21).  Combining  these  equations  yields 

eV-E  +  (Ve)-E  =  p  (3.31) 


Divide  by  e  on  both  sides 

V  •  E  =  —  -  E  =  —  -  (V  In  e)  •  E 
e  e  e 


(3.32) 


Now  take  curl  of  the  third  of  equation  (3.1) 


Vx(VxE) 


(3.33) 


V(V-E)-V2E  =  -  —  (Vx(//H)) 
And  substitute  equation  (3.32)  into  equation  (3.34) 


—  -  (V  In  e)  •  E 
e 


-V2E  =  -  — (Vx(//H)) 


(3.34) 


(3.35) 


One  of  the  laws  for  operations  with  the  div-and-curl  operator  is 

V  x  {(j)F)  =  t/>V  x  F  +  (V  </>)x  F 
and  so  equation  (3.35)  becomes 


(3.36) 


V2E  -  V  —  +  V  [(V  In  e)  •  E]  = 


ju  —  (VxH) 

dr  ’ 


+  -[(V//)xH] 


(3.37) 
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Using  the  fourth  line  in  equation  (3.1),  equation  (3.37)  becomes 


V2E  -  e/j  p-  +  [V(V  In  e)  •  E]  =  V  —  + — [(V//)  x  H] 
C)t  G  Q)t 


(3.38) 


For  nonmagnetic  material,  //  =  juQ .  When  p  =  0  ,  equation  (3.38)  can  be  simplified  to 


<32E 


V2E  +  V  [(V  In  e)  •  E]  -  e/j  — =  0 


dt1 


(3.39) 


Finally,  with  the  help  of  the  refractive  index,  n1  =  ,  and  the  speed  of  light  in  vacuum, 


-2  =  1/ 


c~  =  y_  ,  equation  (3.39)  becomes 


V2E  +  2V[V(lnn)-E]-^r^  =  0  (3.40) 

The  extra  tenn  in  equation  (3.40),  2V  [V(ln  n)  ■  E] ,  will  not  be  zero  when  n  varies  with 

position.  So  electric  fields  Ex,  E  and  E,  may  not  satisfy  the  same  wave  equation. 

Consequently,  we  can’t  make  these  vector  waves  into  a  single  scalar  wave.  Also,  we  will 
get  some  error  by  using  scalar  theory,  even  in  the  homogeneous  medium  when  the 
boundary  conditions  are  imposed  on  a  wave  [20]. 

We  will  consider  scalar  components  and  the  corresponding  scalar  wave  equation 
instead  of  the  full  vector  theory  in  the  thesis. 
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IV.  THE  TRANSFER  FUNCTION 


For  the  purpose  of  our  model,  we  have  to  develop  the  transfer  function  of  a  wall- 
like  medium.  We  assume  that  the  wave  propagates  into  the  medium  with  an  index  of 
refraction  ( nm ).  The  source  is  located  at  r0  =  -z0z,  and  coordinates  are  as  shown  in 
Figure  7. 

The  scalar  wave  equation  with  source  s(rj)  in  the  time  domain  becomes 

V2E-^^  =  s(r,t)  (4.1) 

c  dr 

Using  the  Inverse  Fourier  Transform,  we  can  express  the  source  term  in  terms  of  angular 
frequency  ( co  )  and  oscillation  frequency  ( v  )  as 

s(r,  t)  =  —  [  S(r,  (o)e'Md(o  (4.2) 

In  J-°° 

where  S(r,  co)  is  the  source  tenn  in  the  angular  frequency  domain. 

Using  the  identity  co  =  2nv  and  do  =  2 ndv  ,  equation  (4.2)  becomes 

s(r,t)=  f"  S(r,v)ei2m'dv  (4.3) 

J-  oo 

where  S(r,v)  is  the  source  term  in  the  oscillation  frequency  domain. 

In  the  previous  chapter,  we  showed  that  a  plane  wave  solution  is  in  the  fonn  of  equation 
(3.15).  This  will  satisfy  equation  (4.1),  and  we  can  obtain  the  following  relationship. 


Then  we  can  make  equation  (4.1)  into  a  time  independent  wave  equation  as 

V2E  +  n2k2E  =  S(r,v)  (4.5) 

To  get  the  solution  of  equation  (4.5),  we  put  the  source  at  the  point  r0  =  -z0z  and  use  the 
Green’s  function  which  satisfies  the  following  equation 
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Figure  7.  Point  source  emanating  waves  through  a  medium 


(V2  +  k2n2(z))g(r ,  r0)  =  5{r  +  r0)  =  S(x)S(  y)S(z  +  z0) 


where  k  = 


and  A  is  the  wavelength  in  free-space. 


The  refractive  index  of  our  system  is 


n(z)  =  \ 


n 


1 


for  z  <  0 
for  0<z<L 
for  L<z 


Then  we  can  get  the  solution  for  E  in  the  frequency  domain  as 

E(r,v)  =  \  g(r,r0)S(r,v)d3r0 


(4.6) 


(4.7) 


(4.8) 


Where  g(r,r0)  is  the  Green’s  function  in  the  spatial  domain,  and  the  solution  for  E  in 
the  time  domain  (using  the  Inverse  Fourier  Transform)  as 

E(r,t)  =  I"  E(r,v)ei2xv,dv  (4.9) 

J  -oo 
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In  our  model,  we  assume  that  the  medium  is  essentially  non-conductive,  so  the 
electromagnetic  wave  can  propagate  the  medium  with  little  damping.  Since  we  put  the 

source  at  0  0  along  the  z  direction  we  can  convert  the  equation  (4.6)  into  the 

Fourier  domain. 


G(vx,v  ,z,r0)  =  f  f  g(x,  v,  z,  r0)e  ,i2m'’x+27lVyy)dxdy 

J  -00  J  -CO 

g(x,  y,z,r0)  =  J  G(vx ,  vy ,  z,  r0  y^x+2^d  vdv^ 


(4.10) 

(4.11) 


These  expressions  are  correct  in  our  model,  because  the  medium  remains  constant  in  the 
x  and  v  direction.  Therefore,  there  is  no  changing  of  the  medium  in  these  directions. 

Now  substituting  equation  (4.11)  into  equation  (4.6)  and  using  the  identity  in 
equation  (4.10),  we  obtain 


-(2^)2-(2Wv)2+^y 

oz~ 


G(vx ,vv,z,r,)  +  k2n2(z)G(v  ,vv,z,r,)  =  S(z  +  z0)  (4.12) 


x’  v 


If  we  define 


F  =  2 nvxx  +  2nv  y 


we  get 


|F|==(2w,)!+(2w„)! 

Using  the  wave  number  definition,  we  can  rewrite 


(4.13) 


(4.14) 


k2=k2x+k2y+k2  = 


r  2 

2 

_i_ 

(2^- 

2 

r 2n^ 

9 

“T 

UJ 

i 

It, 

(4.15) 


=  (2xvx)2  +(2xvvf  +(2;rv:)2 


1  1  1 

Vx~T’  Vy~T’  Vz~T 


If  z  ^  — Zq  ,  the  equation  (4.12)  becomes  the  homogeneous  wave  equation  as 
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Figure  8.  Cross  sectional  view  of  the  wave  front 


(4.16) 


(4.17) 


The  solutions  of  the  wave  equations  in  each  region  must  satisfy  the  boundary  conditions, 
which  state  that  the  wave  and  its  derivatives  must  be  continuous  at  all  boundaries.  Then 
we  can  find  the  coefficients  of  each  component  wave  as  they  propagate  through  the 
medium  (Figure  8). 

The  boundary  conditions  are 


gb\ 

=  Gh  1 

dGh 

sgb 

\z=-z0 

\z=-z0 

dz 

dz 

z=~z0 

(4.18) 


Gb\ 

=  Gm\ 

dGb 

_  dG"‘ 

lz=0 

L-=o 

dz 

z= 0 

z=0 


(4.19) 


30 


G"‘ 


\z=L 


(4.20) 


=  GT\  , 

8G'n 

_  8Gt 

\z=L  \z=L 

dz 

dz 

z=L 

We  can  ignore  the  coefficient  AB  of  the  left  going  component  at  z  <  -z0  and  B  r  of  the 
right  going  component  at  z  >  L  ,  because  the  wave  is  outgoing  from  the  source.  In 
addition,  by  imposing  the  Sommerfeld  radiation  condition,  the  surface  integral  at  infinity 
vanishes  [16].  Since  the  waves  are  continuous,  the  Green’s  functions  and  its  derivatives 
also  have  to  be  continuous,  except  at  the  source.  Let’s  find  the  derivative  of  the  Green’s 
function.  We  use  equation  (4.12)  and  integrate  it  on  both  sides. 

J G(yx , vy  ,z,r0)dz  +  f  ( k2n2(z ))  - (2 nvx  f  - (2xvy)2 G(vx , vy  ,z,r0 )dz  =  J S(z  +  z0 )dz 

(4.21) 

The  left  hand  side  becomes  the  first  derivative  of  the  Green’s  function  and  an  integral  of 
the  Green’s  function.  The  right  hand  side  is  the  same  as  the  Heaviside  Step  Function, 
which  is  a  discontinuous  function  whose  value  is  zero  for  negative  argument  and  one  for 
positive  argument.  Then  the  equation  becomes 

Yz  G(vx’vy’z’  ro)  +  J  (k2n2(z)  - (2 7rvx)2  -  ( 27tvy)2)G(vx,vy ,  z, r0)dz  =  H(z  +  z0)  (4.22) 

The  second  term  in  equation  (4.22)  will  vanish,  because  the  Green’s  function  is 
continuous  at  any  given  boundary.  But  the  first  term  will  be  continuous  as  long  as  the 
boundaries  are  z  <  -z0  or  z  >  -z0 .  As  the  first  derivative  crosses  the  source  where  the 
Heaviside  Step  Function  goes  from  0  to  1,  the  difference  across  will  jump  to  1.  Using  the 
same  boundaries  as  the  above,  we  get  the  following: 

If  we  let  the  refractive  index  obey 

nh  (z  <  0) 
n,n  (0  <  z  <  Z) 
nT  ( L  <  z) 


then,  we  can  derive  the  following  equations: 
At  z  =  0 
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Ab  +  Bb  =  A"1  +  B"‘ 


(4.23) 


Ahjk2n;-\Ff  -BbJk2n2b-\F\2  =  Am Jk2n2m-\F\2  -Bm ]k2n2m-\F^  (4.24) 


1  1  ]\Ab]=\ 1  1  ]\A 

X  -XM~X  -Xls 


(4.25) 


where  N,  =  k  n.  -  F~  and  N  =  k  ni  -  F ' 

b  b  mm 


At  z  =  L 


(4.26) 


A2  _l  pP/i*2""  lFl" 


AmJk  nt  -\FYe 


-B"\  k2nl  -\f\  e 


=  AT  Jk2nl  -\F\2 


(4.27) 


-X 


iFFL lr  X 


e1^1  L  0 


(4.28) 


where  NT  =k2n2-\F\~  and 


At  z  =  — zn 


-%B /jzQ'Jk1nl  |F|  jb ^-izoyj^nl  |^|  j^b ^zoyjk2nl-\F\ 


(4.29) 


We  need  to  use  the  jump  condition,  which  is  the  second  one  in  equation  (4.18)  to  obtain 


Ab]k2n2b-\F\ 


2  -izJk2n2b-\Pf 


-Bhjk2n;-\F\ 


2  tJk2n2b-\F\ 2 


=  -BByjk2n2b -\f\ eiz^k2nl~^  +1 


(4.30) 
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o 

T 

i _ 

i 

U‘  1 

’T 

lir 

_ j. 

-jN~beiz°^ 

1 

_ 1 

e~iZoy[Nb 

~BB~ 

"0" 

-jN~beto^b 

0  _ 

+ 

1 

(4.31) 


2  2  |  1 2  •  • 

where  Nh  =  k~nh  -  Fp .  Now  we  can  get  the  solution  for  the  coefficient  of  the  Green’s 
function.  Let’s  write  our  results  in  the  simplified  form 


~Ab~ 

M, 

=  ¥, 

1 

Bh 

2 

~  Am~ 

r 

m3 

B"‘ 

=  m4 

A" 

B" 

A1 

0 


l_ 

_l 

l_ 

1 

~Ab~ 

~BB~ 

"0" 

M5 

=  Mb 

+ 

Bb 

0  _ 

1 

(4.32) 


(4.33) 


(4.34) 


where  the  reflection  and  the  transmission  amplitudes  can  be  expressed  in  terms  of  the 
amplitude  of  the  incident  wave  Ah 


~Ab~ 

Bb\ 

=  m;1m2m~1m4 

1 

^  ° 

1 _ 

'p 

_i 

Ab 

=  /3uAT  Bb=j82 

,AT 

Bb=^Ab  At=  — 

Pi.  i  0u 

Doing  some  algebra,  we  get  the  reflectance  and  the  reflection  coefficient  R  =  r*r . 
where 


(4.35) 

(4.36) 

(4.37) 


Bb  ( 

ftp  +  -  V7i) + (-Rp + v 

^)( 

rN~) 

\ei2L ^ 

r~^~\ 

\[Nt  +  V 

'*■)( 

yj^b  +  \l 

+  \l 

^)( 

\Jnt  —  -\j 

fJT) 

iV  m  ) 
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By  the  same  technique,  we  can  obtain  the  transmittance  and  the  transmission  coefficient 
T  =  t*t. 

where 


_ 4  _ 

(4.39) 


We  also  get  the  other  coefficients  inside  the  medium 


Am 

B'n 


m;‘m4 

p 


at 

0 


Am=P  Ar=^Ab  B"‘  =PltAr  =^Ah 

’  fiu  A,i 


(4.40) 

(4.41) 


A" 


M 

VAX 

vr) 

+  V 

^)(V^+V^)+(-V^+- 

A 

)( 

■\Jnt  -\j 

'A) 

|e«'2£V^r 

2  JNbUNm-JNT 


Bm  _ 

A‘ b  ( Jiw  +  Jn~)  (  Jk  +  Jn~) + (-JnI + Jn~)  (  JnI  -  Jn~)  ei2L^ 

\  v  i  V  7M  7  \  V  b  v  /w  /  \  V  b  V  m  /  \  V  f  V  /w  / 

At  the  source,  we  obtain 


(4.42) 


(4.43) 


5" 

0 


roil 

)[o. 

Uj 

and  if  we  let 


s  =  m; 


i 


^°^(i+v^) 

1 


-fe 


o+VA) 


(4.44) 


(4.45) 


(4.46) 
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then  equation  (4.44)  becomes 


Bb 

0 


5 


(4.47) 


where  BB  =  Xt  lAr  -Su 


A 


Ab 


V  j 


Su  and  0  =  X21AT  -S21 


(4.48) 


We  know  that  the  refractive  indexes  outside  of  the  medium  are  the  same:  Nb  =  NT  =  N . 
What  we  need  for  our  simulation  are  the  amplitudes  of  the  incident  wave  ( Ab  ),  the 
reflected  wave  ( Bh ),  and  the  transmitted  wave  ( Aj  ). 


Ab 


2  4n 


Bb  =r*Ab 


(.N-N.)+(N.-N)e'2Ljii: 


eizo  An 


2 


(4.49) 


(4.50) 


At  =t*Ab 


(4.51) 


Since  the  free  space  Green’s  function  is  a  solution  of  the  Helmholtz  equation  and  satisfies 
the  radiation  condition  at  infinity,  this  should  be  the  solution  in  that  region.  The  following 
are  the  Green’s  functions  we  are  interested  in  : 


For  the  region  -z0  <  z  <  0 , 


G(vx,vy,z,rn)  =  < 


eiz^  bfk 

- -e  y 


2y[N 

<N-NJ+(N:-N)ei2L^ 


right  going 

Jz0  -Jn 


N  +yJN)  -  JN,- 


N)2eaL ^ 


e  '-A2  lrl  left  going 


(4.52) 


For  the  region  z  >  L  , 
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G(vx,vy,z,r0)  = 


h4n 


N+JN  I  e^N"'  -UN.  -ylN  I  e 


iLjN, 


ei:°^  ■ 

e'-v  II  right  going 


2  4n 


(4.53) 


We  can  also  get  the  coefficient  Am  and  B"“  by  solving  this 


Am 

Bm 


M~lMx 


Ah 

Bh 


(4.54) 


For  the  region  0  <  z  <  L  ,  we  get  these  coefficients 


Bm 


!Z0VV 


2\ 

fN  | 

K- 

-Jn) 

\ei2L V 

'aT 

[4n +^j 

2 

-1 

(v 

4n] 

|2 

Finally,  we  obtain  the  Green’s  function  in  the  media 


G(vx,vy,z,r0)  =  A'ne‘ 


right  going  wave 


G(vx,y  ,z,r0)  =  B,ne  ^  ^  left  going  wave 


(4.55) 


(4.56) 


(4.57) 

(4.58) 
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V.  THE  ANGULAR  SPECTRUM 


A.  TRANSFORM  SPHERICAL  WAVES 

The  angular  spectrum  method  is  a  technique  for  modeling  the  propagation  of  a 
wave  field.  This  method  involves  expanding  a  complex  wave  field  into  a  summation  of 
many  plane  waves.  As  the  spherical  waves  propagate  through  a  wall,  they  are 
appropriately  described  by  a  planer  expansion,  because  waves  propagate  in  all  directions 
spherically,  and  the  boundary  is  plane.  We  have  to  transfonn  the  spherical  waves  to  plane 
waves,  since  the  reflected  wave  from  a  plane  media  with  a  spherical  wave  is  also 
spherical  [21].  To  obtain  the  angular  spectrum  g(x,y,z,r0 )  we  need  to  take  the  Inverse 

Fourier  Transform  of  the  right  going  wave  in  the  region  -z0  <  z  <  0  in  equation  (4.52). 
We  obtain 


g(x,y,z,r0)  =  J  J  G(vx,v  z,r0)e‘(2m'lX+2^y>dvxdv 


(5.1) 


In  this  region,  the  refractive  index  is  n  =  1 ,  and  equation  (5.1)  becomes 

izo  yjk2  “(2  Wv  )2  -(2  XV y  )2 


poo  poo 

g(x,y,z,r0)  =  J  J 


e 

o  Ik 2 


/"A2  (2 jzv  )2  -( 2ttv  )2  i2x(v  x+v  y)  .  , 

2  v  e  ”  dvdv. 


2\lk  ~(2xvx)  ~(27rvy) 


(5.2) 


ik(z(s+z)di-(2nvxlk)2-(2jrvylk)2 

g(x,y,z,r0)  =  J  ^  J  - 1  =e{2m,’‘x+2m,^dvxdvv  (5.3) 


2kx\\ -  (2 nvx  /  k)1  -  (2 nvy  /  k) 


—  2n 

Let’s  suppose  a  wave  vector  k  has  magnitude  — j-  and  direction  cosines  (a,  /3,y)  ,  (See 
Figure  9)  then  k  is  expressed  as 


k  =  kxx  +  k  y  +  k.z  =  k(ax  +  J3y  +  yz) 


2  n 

T 


(ax  +  fdy  +  ^1  -a2  - /32 z) 


(5.4) 

(5.5) 
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since  direction  cosines  are  a'  +  p2  +  yn  =1  where 


kx  =  k  ■  x  =  ka 
ky=k-y  =  kj3 
k,  =  k  ■  z  =  ky 


Then  we  can  get  the  following  identities 


kx  =  ka,  k  =k/3  =>  a  = 


2  nvr 


2  nvv 

fi  =  ~rL 


da  = 


2_n£_  _  dpJ_^L 


Using  these  identities,  equation  (5.3)  becomes 


g(a,p,z,r(])  =  \  J 


OO  <•  CO 

-00 


ik(zt]+z)-J\  a2  (f~ 


ik(ax+/3y) 


8  7T2 


2kyj\-a2-p2 

\  .00  1-00  Jky(z0+z)  ik(ax+Py) 

f  f  — 

I  J  -00  J  -00 


f  k ) 

da 

f  k  3 

[lx) 

\2n  j 

dp 


Y 


-dad  P 


From  H.  Weyl’s  expansion  of  the  spherical  waves  into  planar  waves  [22] 


Figure  9. 


The  wave  vector  k 


(5.6) 
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(5.7) 


iK\r—r  .  -j 

f - =  _L  r  ['  1 -^°i*-M(y-y'W*-?)])dad P 

\r  —  r'\  Ik  ^_0°  ^~°°  y 
where  r  =  xx  +  vv  +  zz,  r'  =  if  +  y'y  +  z'z 
Substituting  equation  (5.6)  into  equation  (5.7),  we  get 


g(x,y,z,r0)  = 


7  1 

k  e 


i\n  \r-rn 


and  applying  the  far- field  approximation  to  equation  (5.8),  we  obtain 

keikr 


g(x,y,z,r0)n- — e~‘  'r° 
lAnr 


(5.8) 


(5.9) 


B.  REFLECTIVITY  OF  MEDIA 

In  the  previous  chapter,  we  obtained  planer  expansions  of  Green’s  functions.  To 
calculate  the  reflectivity  on  the  surface,  we  need  to  employ  spherical  coordinates.  In 
equation  (4.52),  we  see  the  reflectivity  ( R  )  from  left  going  waves  as 

*=  (N-N.MN.-N)e^  (5i0) 

(skL+'tu)  -(sIK.-'ki) 

where, 

4Tm=jnik2-Q.KVx)2-Q.KVyf  ,  4n  =  ^Jk2  -(2kvx)2  ~(2kvv)2  (5.11) 

The  spherical  expansion  of  the  reflected  wave  can  be  obtained  by  converting  R  into 
spherical  coordinates.  Using  the  relationship  in  equation  (5.4)  and  the  spherical 
coordinate  representation,  we  obtain 

2 kvx  =  k sin 6 cos ^ ,  2 kv  =ksmOcos<j),  2 nv_  =  k cos 6  (5.12) 

Substituting  equation  (5.12)  into  equation  (5.1 1),  we  get 
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(5.13) 


■sjNm  =  kyjnl  -  sin2  0 , 4n  =  W  1-sin2  0  =  k  cos  6 


And  inserting  these  results  into  equation  (5.10),  we  get  the  reflectivity  as 


R{0)  = 


(}-n2m)  +  (n2m-l)e 


2kL\fn 


( \]nl  -  sin2  6  +  Vl  -  sin2  9  j  -sin2  <9  -V  1-sin2  6 

/I  2  \  .  (  2  i  \  i2kLJnlj  -sin2  6 

_ (t-flJ  +  K-  l)e  v  _ 

■sjnl  -  sin2  0  +  cos^  j  - -  sin2  6  -  cos  0  \  e'2tt'/""“sin  6 


i2kLyji 


izkL\  n„-  sin 


(5.14) 


(5.15) 


C.  EVANESCENT  WAVES 


The  electric  field  does  not  abruptly  vanish  at  the  boundary,  because  the  phase 
difference  between  the  incident  and  reflected  waves  prevents  the  complete  destructive 
interference  required  to  eliminate  the  transmitted  wave.  The  transmitted  field  that  extends 
beyond  the  boundary  of  the  dielectric,  when  a  wave  is  totally  internally  reflected,  is 
known  as  the  evanescent  wave  [19]. 

If  we  choose  a  coordinate  system  in  which  the  boundary  lies  in  the  x-y  plane 
and  k  lies  in  the  x  —  z  plane  and  suppose  we  have  a  transmitted  wave  of  the  form 


E,  =  E„e 


i{kt-r-cot) 


then,  we  can  rewrite  this  with  Snell’s  law  as 


kt-r  =  ktx  sin  Qt+ktz  cos  0t 


=  k, 


xsin#  .  sin#  , 
- l  +  iza  — ^--1 


(5.16) 


(5.17) 

(5.18) 


Substituting  equation  (5.18)  into  equation  (5.16),  we  get  the  transmitted  field 

Et  (r,  t )  =  El()e-k'zel(k'x-at)  (5.19) 
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where  r  is  a  constant  and  k'  =  kt  sin  0t .  This  wave  propagates  along  the  boundary  with 
propagation  constant  k' ,  diminishing  exponentially  with  constant  r . 


Let’s  look  at  the  region  -z0  <  z  <L  and  examine  how  much  the  wave  evanesces. 
We  need  equation  (4.52)  and  equation  (4.57)  for  the  reflection  and  transmission. 


g(x,y,z,r0)  =  A(vx,v 


V  -i(z+Z0  )yjk2  -(2 7TVX  )2  -(2 7TVy  )2  i(27TVxX+27TVyy) 

,  )C  c 


(5.20) 


g(x,  y,  z,  r0)  =  i'(vx,vj/z+Zo)^2‘(2^)2‘(2^)2e,'(2^+2^) 


(5.21) 


For  the  reflected  wave  in  equation  (5.20),  there  is  an  angle  of  incidence  that  results  in  a 
transmission  angle  that  is  parallel  to  the  surface.  If  the  incident  angle  increases  over  the 
critical  angle,  the  transmitted  wave  will  evanesce.  In  the  region  z  <  0 ,  the  intensity  of 
reflection  wave  is 


I  =  g  g  =  A  Ae 


2(z+zt))^(27cvxx)1-(27ivyy)2  - k 2 


(5.22) 


As  z  goes  negative,  the  intensity  becomes  smaller,  so  the  evanescing  wave  will  be  much 
smaller.  Therefore  in  this  model  system,  we  will  ignore  the  evanescent  waves  as 
unimportant. 


n 


Figure  10. 


The  wave  vector  kt,  kr  and  kt 
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Figure  1 1 .  Reflection  (polystyrene,  quartz,  glass) 


D.  MATERIAL  ANALYSIS 


Typical  materials  of  common  walls  are  wood,  tile,  and  concrete.  These  materials 
have  their  own  refractivity  defined  by  the  index  of  refraction 


c  sin  6, 

n  =  —  = - - 

v  sin  0r 


(5.23) 


The  refractive  indexes  of  materials  depend  on  the  frequency  with  which  it  is  measured. 
Unfortunately,  indices  of  refraction  in  the  terahertz  region  have  not  been  well  tabulated. 

Not  all  light  striking  a  transparent  material  is  refracted  as  shown  in  Figure  10  [23].  The 
reflectance,  R  is  related  to  the  index  of  refraction  by 
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R  = 


lk+ 1 


V'b+i 


\2 

~nk 

+  fh; 


(5.24) 


(when  the  incident  angle  6f  =  0 ). 

From  equation  (5.24),  we  can  predict  that  the  higher- n  materials  are  the  more  refractive. 
In  this  section,  we  will  analyze  materials  using  MATLAB  simulations  based  on  equation 
(5.15)  and  find  appropriate  frequencies  that  can  penetrate  wall  materials  without  much 
loss  and  with  good  reflection  of  objects.  For  this  simulation,  we  assume  that  we  are  using 
plane  waves  and  nonnal  incidence. 

One  of  the  most  attractive  advantages  about  terahertz  frequencies  is  their  ability  to 
penetrate  many  common  materials.  But  the  terahertz  pulse  can’t  penetrate  materials  with 
high  water  content,  so  we  assume  the  materials  used  in  this  simulation  are  pre-dried. 
Table  II  shows  index  of  refraction  for  common  materials  in  terahertz  frequency  [4,  8]. 

Figure  1 1  shows  that  as  we  increase  the  index  of  refraction,  the  reflection 
becomes  large.  In  Figure  12,  we  used  refractive  indexes  of  common  wall  materials.  In  the 
case  of  tile,  rock,  and  sand  the  average  reflection  is  quite  high.  But  if  we  look  at  this 
figure  closely,  we  can  find  the  frequencies  that  have  low  reflection  from  these  materials. 
These  frequencies  are  at  about  0.94THz,  1.44THz,  etc.  and  can  be  used  for  Through-the- 
Wall  Imaging. 
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Reflection 


Polystyrene  foam 

1.01 

Ochroma  pyramidale 

1.08 

Lophira  lata  (wood) 

1.49 

Skin 

1.4 

Body  fluid 

3.9 

Sand 

1.67 

Glass 

1.8 

Quartz 

1.5 

Tile 

1.6 

Stone 

1.6 

Table  2.  Index  of  refraction  for  common  materials. 
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Figure  12.  Reflection  (lophira,  tile,  sand) 
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We  also  used  the  index  of  refraction  of  human  skin  and  body  fluid.  Figure  13 
illustrates  the  reflection  from  body  fluid  is  very  high  at  frequency  0.94THz.  Consequently, 
if  we  choose  0.94THz  frequency  for  Through- the- Wall  Imaging,  we  will  get  more  data 
from  the  human  body  and  get  a  good  image. 
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Figure  13.  Reflection  (body  fluid,  skin) 
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VI.  INTEGRAL  EQUATION 


A.  WAVE  PROPAGATION  EQUATION 

In  this  thesis  we  are  assuming  that  the  complex  dielectric  constant  and  magnetic 
penneability  vary  slowly  over  a  wavelength  of  the  electromagnetic  wave,  and  we  are 
neglecting  polarization  effects.  Then  we  can  write  the  Helmholtz  equation  as 

V2£(r)  +  £24(r)£(r)  =  0  (6.1) 

If  we  add  k2E(r)  on  both  sides  we  get 

V2E(r)  +  k2E(r)  =  -k2(nobj(r)2  -l)E(r)  (6.2) 

0(  r) 

where  0( r)  is  the  target’s  object  function  and  is  denoted  by 

0(r)  =  k2(n2ibj-\).  (6.3) 

Then  the  total  electric  field  E{ r)  solution  to  equation  (6.2)  is 

E{  r)  =  Einc  (r)  +  £g(r  -  r0)O(r0)£(r0  )dr0  (6.4) 

where  g  is  the  Green’s  function  of  the  Helmholtz  equation  and  is  the  solution  of  the 
differential  equation 

(V2  +  k2)g(r  - r0)  =  -S(r  -  r#)  (6.5) 

where  £(r-r0)  is  the  Dirac  delta  function.  In  our  model  system,  the  Green’s  function 

varies  as  the  observation  point  changes.  From  the  previous  chapter,  we  found  Green’s 
functions  in  the  spatial  frequency  domain  as  Eq.(4.52),  Eq.(4.53),  Eq.(4.57),  and 
Eq.(4.58).  So  in  order  to  get  g(r-r0),  we  have  to  transform  these  Green’s  functions 
using  the  Inverse  Fourier  Transform  as 

g(r  ~ro)  =  JJJ  G(v x,v y,z,ra)e,2*<v'x+v‘yi  dv  xdv y  (6.6) 
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We  denote  each  Green’s  function  by  their  domain  as 

&/(r-r0)  for  z<  0 
&//(r~ro)  for  0<z<L 
S-///(r  -ro)  f°r  L<z 

and  we  will  use  g/H(r-r0)  for  our  model  system  in  the  following  section. 

B.  THE  FIRST  BORN  APPROXIMATION 

The  Born  approximation  is  perhaps  the  simplest  and  most  widely  used  EM 
scattering  approximation.  It  is  based  on  the  assumption  that  scattered  electric  fields  inside 


48 


scatterers  are  negligible  compared  to  the  normal/background  electric  fields.  Hence  the 
second  term  on  the  right  hand  side  of  equation  (6.4)  is  small  in  comparison  to  the  first 
tenn  within  the  integral.  In  the  integrand  can  approximate 

E(r)*Einc(r )  (6.7) 

Therefore  the  total  electric  field  in  equation  (6.4)  is  calculated  as 

E(r)  =  E,nc(r)  +  \y  g(  r  -  r0)O(r0)Einc(r0)d\  (6.8) 

Consequently,  the  scattered  field,  Es ,  can  be  identified  with 

Es  =  \y  g(r  -  r0)O(r0)Einc  (r0)d\  (6.9) 

C.  THE  FIRST  RYTOV  APPROXIMATION 

The  Rytov  approximation  assumes  that  the  incident  wave  perturbation  caused  by 
the  target  can  be  described  by  a  change  of  phase  in  the  reference  wave.  The  Rytov 
approximation  is  obtained  by  representing  the  total  field  as  a  complex  phase  [22] 

E(  r)  =  em  (6.10) 


We  consider  the  scattering  of  a  scalar  wave  in  an  inhomogeneous  medium.  The  field 
satisfies  the  reduced  wave  equation 


V2E(r)  +  k2E(r)  =  0 

(6.11) 

Substituting  equation  (6.10)  into  equation  (6.1 1),  we  get 

vV(r)+cV(r)  =0 

(6.12) 

Using  the  identity 

yV(r)  =V2</>em +(V</>)(V</>)em 

(6.13) 

we  substitute  equation  (6.13)  into  equation  (6.12)  and  divide  by  e*(r) 

get  the  following  Riccati  equation 

on  both  sides  and 

V2</>  +  (V</>)(V</>)  +  k2  =  0 

(6.14) 

We  denote  the  refractive  index  as 

k(r)  =  k0n(  r)  =  k0(  1  +  nS(r)) 

(6.15) 
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The  zero  perturbation  field  for  n(r)  =  1  can  be  written  as  E.( r)  =  Ulr)  ,  and  its  phase  will 
satisfy  equation  (6. 14) 

VV,+(V$)2+*o2  =0  (6.16) 

We  can  express  the  total  complex  phase,  (j) ,  as  the  sum  of  the  incident  phase  </>.  and  the 
scattered  complex  phase  <f>s 


*  =  *+*,  (6- IV) 

Now  substitute  equation  (6.17)  into  equation  (6.14) 

V2(4+4)  +  (V(4+4))2+*2  =  0  (6.18) 

Expanding  this  equation 

V2^.  +  V20s  +  (V^)2  +  (V^)2  +  2(V^.)(V^)  +  k2  =  0  (6.19) 

and  subtract  equation  (6.16)  from  equation  (6.19),  we  get 

V2^  +2 .(V4XV*,)  +  (*2  -k2)  =  -(V</>s)2  (6.20) 

V2^.s  +  2(V <pt )(V <ps)  =  -((V <t>s )2  +  kl (n2  - 1))  (6.21) 

Using  the  identity 

V2(£,4)  =  <i>sV2Ei+2(VEi)(V<i>s)  +  EiV2<i>s  (6.22) 

and  assuming  a  plane  wave  for  the  incident  field,  so  that 

E.  =  Aeik°sr  =  Ae¥l(r)  (6.23) 

then  we  find  the  following  expression 

V2Ei=-k2Ei  and  VE  =  ES/fr  (6.24) 

So  equation  (6.22)  becomes 

2  EFfi -V^+E^  =V\E^a)  +  k^s  (6.25) 

Now  substituting  equation  (6.25)  into  equation  (6.21)  ,  we  get 

V2(£,4)  +  *2(£,4)  =  -E,  {(V4)J  +kl(n2  -1)}  (6.26) 

Solving  this  equation  using  the  free  space  Green’s  function,  we  get 

4(r)4(r)  =  j  S(r  -r,)£,(r,){(V4(r,)2  +  *2(n>,)  - 1)}*,  (6.27) 


'Ain-  '  f.t'(  r  -  r„  )£('■„  IC  -  '  f.t-lr-r,)i'.(,-ll)/vl«J(,,)-l)rfr11  (6.28) 

E,(r)J  E,(r)J 
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By  the  Rytov  approximation,  the  first  tenn  in  equation  (6.28)  is  very  small,  so  equation 
(6.28)  can  be  approximated  as 

«t-(r)  =  V7Tl«<r“r«)£.(r»)*:»<',I(,'»)-1)‘,r»  (6'29) 

E,{  r)J 

If  we  denote  the  target  object  function  0(r )  as 

0{r )  =  kl  (n2  (r0 )  - 1)  (6.30) 

Then  equation  (6.29)  becomes 

&(*0  =  T7T  f  g(r  -r0)£,(r#)O(r0)Jr0  (6.31) 

E,(  r)J 

The  corresponding  solution  for  the  total  field  is 

E(r)  =  Et(r  )e*>  (6.32) 

But  it  is  only  valid  when  (V^s)2  □  £2(«2(r0)  -1) 

From  equation  (6.10),  the  total  electric  field  is  given 

E(r)  =  Ei(r)  +  Es(r)  =  eMr)+Mr)  (6.33) 

Rearranging  this  equation  for  the  scattered  field,  Es  (r)  ,  we  get 

Es(r)  =  E. (r )(e^ (r)  - 1)  (6.34) 


Substituting  equation  (6.31)  into  equation  (6.34),  finally  we  get  the  scattered  field  as 


Es(r)  =  Ei(r0) 


— *-r  f  g(r-r0  )E,  (r„  )O(r0  >/r0 

eE‘(r) 


-1 


(6.35) 


Note  that  when  the  argument  of  the  exponential  is  small,  then  we  can  approximate 

'  Tjs<r-r0)£-,(r0)O(r„)fl'rll  ^  } 


Mr) 1 


1  +  7T7  f  “ r« )Et ^)O(r0 )dr0 

E,(  r)J 


(6.36) 


and,  in  this  case,  the  Rytov  approximation  reduces  to  the  first  Born  approximation. 
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VII.  SIMULATION 


A.  RECONSTRUCTION  ALGORITHM 


For  many  reasons,  it  is  very  useful  to  convert  functions  in  the  time  domain  into 
functions  in  the  frequency  domain.  The  primary  reason  is  that  it  gives  us  an  easy  way  to 
solve  difficult  problems.  For  instance,  when  we  meet  the  convolution  integral  in  the  time 
domain,  we  can  express  this  as  a  simple  multiplication  integral  in  the  frequency  domain. 

To  apply  the  Rytov  approximation  for  object  reconstruction,  we  need  to  estimate 
the  scattered  field  Es( r),  which  is  denoted  in  equation  (6.35).  From  equation  (6.35)  we 
have 


Zf(r)ln 


£,(r) 

£,(r o) 


g(r-r0)£,(r#)O(r0>/r0 


(7.1) 


Since  a  kernel  is  formed  as  a  convolution  between  the  Green’s  function,  the  scattering 
potential,  and  the  incident  electric  field,  we  make  a  reconstruction  algorithm  in  the 
frequency  domain  by  Fourier  Transform  of  the  measurement  data.  The  left  term  in 
equation  (7.1)  is  our  measurement  data.  If  we  denote  this  by  m(r ) ,  and  then 

m(r)  =  f  g(r  -r0)£,(r0)O(r#)dr#  (7.2) 

J  VR 

We  extend  the  volume  integral  to  infinity  and  take  2D  Fourier  Transform  by  assuming 
receivers  will  be  place  at  fixed  points  (xi,yi,z).  Then  equation  (7.2)  becomes 


where, 


M(v„vv)  =  G(v  ,v  )X(v  ,v  ) 


x(y*’vy)  =  \  J  Ei (xo , To )O(x0 , y0 )e  ‘lKiv*x+Vfy}dxdy 


(7.3) 


(7.4) 


By  using  equation  (7.3)  and  applying  equation  (2.9)  we  can  get  X(vx,vy)  as 
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=  [(G{vx,vy)TG{  vx,vyj)  1  G(vx,vyY  j  M(i/r,i/v) 


(7.5) 


Substituting  equation  (7.5)  into  equation  (7.4),  we  can  get  the  “best  object”  O(r0 )  and 
taking  the  Inverse  Fourier  Transfonn  of  this  as 


0(r°)  =  yr^i  I  \  \(G(V*>Vyy  G(V*>Vy))  1  G(Vx,Vy)T\M(Vx,Vy) 
^ i  vro  ' 


i27i{vxx+v  y) 


dxdy  (7.6) 


Finally,  we  get  the  “best  object”  O(r0 )  in  the  spatial  domain.  The  incident  wave 
Ei(r)  =  Ae~,kz  becomes  E'(r)  =  tAe  lkz  after  passing  through  a  transparent  media  with 
transmittance  t .  The  Green’s  function  to  be  used  in  equation  (7.6)  is 


G(vx,vy,z,r0)  = 


4^N^te-iL^ 


N  +  JN, 


-iLJNm 


e'0^  die2- |if 

2 


(7.7) 


B.  MATLAB  CODE  IMPLEMENTATION 

The  total  field  at  the  receivers  will  be  composed  of  the  incident  field  from 
transmitters,  the  reflected  field  from  the  wall,  the  scattered  field  from  the  target,  and  the 
evanescing  field.  In  the  previous  chapter,  we  showed  the  evanescent  field  is  negligible  in 
comparison  with  the  incident  field. 

For  our  image  processing,  we  are  only  interested  in  the  scattered  field  data  and 
use  equation  (6.35)  and  equation  (7.6)  to  create  scattered  field  images.  Due  to  the 
limitation  of  computer  ability,  we  are  restricted  in  the  number  of  scattering  points  to  the 
sampling  frequency,  and  bandwidth  to  be  examined. 
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(b) 

Figure  15.  Amplitude  of  field  scattered  from  5  points  using  Born  (a)  and  Rytov  (b) 

approximations 

Figure  15  shows  5  scattering  points  applying  the  Born  and  Rytov  approximations 
at  z  =  2m  .  We  will  add  Gaussian  noise  to  this  data  and  use  Tikhonov  regularization, 
truncated  SVD  methods,  and  the  L-curve  method  discussed  in  chapter  2.  All  MATLAB 
codes  are  presented  in  Appendix  D. 
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c. 


SIMULATION  RESULTS 


To  compare  the  Born  and  Rytov  approximations,  first  we  simulate  with  equation 
(6.9)  and  equation  (6.35)  separately.  We  get  Figure  15  (a)  from  the  Born  approximation 
and  Figure  15  (b)  from  the  Rytov  approximation.  These  two  figures  look  very  similar, 
and  they  localize  the  scatterers  equally  well.  But  the  strength  of  each  peak  differs  by 
1.2%  in  both  figures.  The  difference  is  formed  by  subtracting  equation  (6.35)  from 
equation  (6.9)  and  simulating  again.  Figure  16  shows  us  the  differences  between  them. 
The  differences  of  each  peak  are  approximately  0.16%,  and  both  approximations  perform 
equally  well. 

D.  SIGNAL  TO  NOISE  RATIO 

For  the  simulation  of  measuring  a  scattered  electric  field,  we  have  to  add  noise  to 
the  original  signal  from  the  Rytov  approximation,  since  signals  are  corrupted  by  noise  in 
nature.  We  will  use  the  Gaussian  noise  that  is  commonly  used  as  an  additive  noise  model 
in  physics.  Since  the  Signal  to  Noise  Ratio  (SNR)  is  defined  as  the  ratio  of  signal  power 
to  noise  power,  we  calculate  the  SNR  by  taking  the  ratio  of  the  sum  of  the  signal  and  the 
Gaussian  noise  in  the  frequency  domain.  Furthermore  we  find  the  minimum  SNR  to 
recognize  the  image  of  objects  by  increasing  a  noise  factor. 
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Figure  16.  Differences  between  Born  and  Rytov  approximations 


Noisy  Data  with  SNR=  -20.0  dB 


Figure  17.  Scattering  points  with  Gaussian  noise  SNR=-20dB 
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E. 


TIKHONOV  REGULARIZATION 


Figure  17  shows  the  scattered  signals  with  noise.  We  use  Tikhonov  regularization 
to  mitigate  these  noise  effects.  During  this  process  there  is  a  trade-off  between  the 
regularized  output  and  the  original  data.  In  order  to  optimize  the  trade-off,  we  have  to 
choose  the  best  regularization  parameters.  To  find  the  regularization  parameters  we  will 
use  L-curve  method.  The  L-curve  method  is  based  on  a  log-log  plot  of  the  nonn  of  the 
residual  versus  the  solution  for  a  range  of  values  of  regularization  parameters  [14].  The 
common  expression  of  the  Tikhonov  regularization  is 

Fa=\\AX-b^2+a\xf2  (7.8) 


II  m2  .  I.  .|2 

where  ||AU-b||]  is  a  residual  norm,  \\X\\  is  a  regularized  solution,  and  a  is  the 

regularization  parameter.  The  “horizontal”  part  of  the  L-curve  is  dominated  by 
regularization  errors  that  occur  from  over-smoothing  and  the  “vertical”  part  is  dominated 
by  perturbation  errors  that  occur  from  under-smoothing.  Hence,  we  choose  our  optimal 
solution  as  the  one  at  the  corner  of  the  L-curve  in  Figure  18(b).  So  a  particular 
regularization  parameter  is  a  =  4.6 .  Using  this  parameter  we  can  get  a  good  recovered 
image  in  Figure  19.  However,  Figure  20  shows  us  that  if  we  use  regularization  parameter 
smaller  than  the  optimal  one,  we  will  get  an  unstable  image. 


120  - 


100  - 


(a) 


alpha  VS  residual  norm 


X:  4.614 
Y:  36.22 


O' — 1 - L 


6 

alpha 


f - r- 


10  12 
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L-curve  for  Tikhonov  Regularization 


Figure  18.  Regularization  parameter  (a)  and  L-curve  (b) 
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F. 


TRUNCATED  SVD 


Since  the  operator  A  of  TWI  has  large  singular  values,  we  can  apply  the 
truncated  SVD  method.  Using  the  singular  value  spectrum  plot  in  Figure  21,  we  can  find 
values  /C  =  0.34  and  ,/C  =  0.10  where  the  singular  values  change  rapidly,  and  truncate 
with  these  values.  Figure  22  is  obtained  by  truncating  /C  =  0.34  ,and  Figure  23  is 
obtained  by  truncating  A3  =  0.10  .  These  two  figures  are  good  examples  of  optimal 
truncation  parameters. 

On  the  contrary,  Figure  24,  obtained  at  JC  =  100  ,  shows  that  a  too-large  JC  affects 
the  accuracy  of  estimation. 


Singular  Values  Spectrum 


Figure  2 1 .  Singular  value  spectrum  of  transfer  operator 
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-2 


TSVD  w/  truncation  at  kappa=  0.34 
SNR=-20.0  dB 


Figure  22.  TSVD  K  =0.34  at  SNR=-20dB 
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TSVD  w/  truncation  at  kappa=  0.10 
SNR=-20.0  dB 


Figure  23. 


TSVD  K  =0.10  at  SNR=-20dB 


TSVD  w/  truncation  at  kappa=  100.00 
SNR=-20.0  dB 
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Figure  24.  TSVD  K  =100  at  SNR=-20dB 
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VIII.  CONCLUSION 


In  recent  years,  a  variety  of  techniques  for  THz  imaging  have  been  developed  for 
the  detection  of  weapons,  explosives,  and  drugs  at  security  check  points  and  for  the 
detection  of  illegal  immigrants  at  the  border.  This  thesis  is  a  review  of  some  of  this  work. 
We  think  our  model  may  be  useful  for  military  and  security  officers  to  detect  enemies 
hiding  in  a  room  or  behind  a  wall. 

In  this  thesis  we  have  applied  the  Bom  and  Rytov  approximations  to  the  problem 
of  imaging  from  scattered  field  measurements.  Also  we  examined  the  Tikhonov 
regularization  parameter  selection  with  the  L-curve  method.  The  simulation  result 
obtained  by  the  Rytov  approximation  has  been  compared  to  one  yielded  by  the  Bom 
approximation.  The  difference  between  them  is  insignificant.  Therefore  the  Rytov 
approximation  appears  to  be  no  better  than  the  Born  approximation.  It  would  be  useful  to 
apply  this  simulation  to  real  data  in  the  future.  In  general,  the  Rytov  approximation  is 
more  accurate  than  the  Bom  approximation  especially,  at  higher  frequencies  [24]. 

Further  modeling  efforts  are  also  required.  We  have  assumed  the  medium  (wall) 
is  linear,  nonmagnetic,  and  nonconductive.  But  in  reality,  wall  materials  are  not  exactly 
like  this.  Some  materials  are  nonlinear,  magnetic,  and  conductive.  So  for  this  case,  the 
wave  equation  has  to  be  adjusted  to  derive  a  new  Green’s  function  to  fit  this  model,  and 
we  have  to  apply  vector  theory  instead  of  scalar  theory  to  the  wave  equation. 

For  more  accurate  image  reconstmction,  we  also  need  to  consider  dispersive 
effects  of  different  kinds  of  walls  and  have  to  include  multiple  scattering  effects  to  the 
Lippman-Schwinger  equation. 
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APPENDIX  A.  SINGULAR  VALUE  DECOMPOSITION 


Let  KeO  MxN  is  a  linear  operator  which  maps  vectors  /  e  v  to  vectors  d  e □  M  . 
We  will  consider  two  symmetric  matrices  KrK(jV  x  jV)  and  KKr(MxM).  Since  these 
matrices  are  square  and  symmetric,  we  can  get  eigenvectors  and  eigenvalues  of  each.  If 
/Lf  is  the  eigenvalue  and  f.  is  the  eigenvector  of  KrK  . 

The  eigen  equation  for  KrK  is 

K'K/: (A.i) 
Applying  the  same  method,  we  get  the  eigen  equation  for  KKr 

KKrd,.  =  yidi  (A. 2) 

If  we  take  the  eigenvector  f.  of  KrK ,  and  the  eigenvalue  then  we  get  vector 

K/l  5*  0  ,  this  can  be  proved  here 

(KKr)K/i  =K(KrK)/i  =  K(Ai)fi  =Aj(Kfj)  (A.3) 

Hence  K/i  is  an  eigenvector  of  KKr  and  each  eingenvalue  ^  0  of  KrK  must  be  an 
eigenvalue  of  KKr  .  Similarly,  this  is  also  true  each  eigenvalue  yt  ^  0  of  KKr  must  be 
an  eigenvalue  of  KrK  . 

If  there  are  r  non-zero  eigenvalues  of  KrK  (or  KKr),  where  r  <  m  and  r  <n ,  we  can 
express  it  like  this 

\=7i 
/l2  =  y2 

K  =  Yr 

A+l  =  0  Yr+1  =  0 

K  =  0  Ym  =  0 
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Therefore  the  linear  equation  is  expressed 


/  =  ¥Jd  d  =  Kf 


(A. 4) 


Then  we  normalize  eigenvectors 


/  = 


K  Td 
YJd\ 


The  singular  values  <rk  of  K  are  defined  by 

Wk  I  =  ||KX  l  =  (jk  =  \[\=4n 

So  equation  (A. 5)  becomes  for  K<M,N. 

KX  =  °kfk  »  K/,  =  <7kdk 

but  for  K  >  M,  N ,  the  equation  becomes 

Kfk  =  0  for  k  =  r  +  1,--- N 
¥fdk=0  for  k  =  r  + 


(A. 5) 

(A. 6) 

(A. 7) 

(A. 8) 
(A. 9) 


On  equation  (A. 7),  we  take  K f  =  <jkdk  and  multiply  by  f'  on  both  sides.  Then  it 
becomes 


K  r=i,aM 

i= 1 

And  the  other  one  in  equation  (A.7)  becomes 

K  =  Yj(T,difT 
1=1 

We  can  write  equation  (A.  1 1)  in  matrix  form  as 


0 

0" 

-...  fx  ... 

K  = 

^2 

0 

o-2 

0 

...  f2  ... 

0 

0 

(A.  10) 


(A.  11) 


(A.  12) 


or  K  =  DUFt 

where,  U  is  a  diagonal  matrix  of  the  singular  values,  and  D  and  F  consist  of  the 
eigenvectors  corresponding  to  those  singular  values.  /  and  d  are  the  orthonomal  vectors 
and  called  right  singular  vectors  and  left  singular  vectors  respectively. 
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We  can  make  the  inverse  of  K  as 

K  1  =  (DUFT)~l  (A.  13) 

We  know  FT  and  D  are  orthogonal  matrices,  because  they  are  composed  of 
eigenvectors  which  are  orthogonal.  Equation  (A.  13)  becomes 

K  1  =  {FT)~lU~lDl  (A.  14) 

K-1  =  FU~xDt  (A.  15) 

U~l  is  a  diagonal  matrix  composed  of  — .  Since  singular  values  of  K  are  in  the  order  of 

<j 

<7,  >  <r2  >•••>  <7r  >  0  (A.  16) 

elements  in  a  matrix  U  1  will  increase  infinitely  [10].  Hence,  the  norm  of  inverse 
mapping  K  1  also  will  increase.  This  is  why  we  need  a  regularization  method  for  the 
inverse  imaging. 
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APPENDIX  B.  GREEN’S  FUNCTION 


Green’s  function  satisfies  a  wave  equation  driven  by  a  point  source. 


(V2  +k2)i//(r)  =  S(r) 

(B.l) 

To  obtain  a  solution  to  wave  equation  need  to  seek  the  Green’s 

solution  to  the  following  equation. 

function  that  is  the 

(V2  +  k2)G(r,r')  =  —S(r  —  r') 

(B.2) 

If  we  let  arbitrary  source  S ( r ) 

S(r)  =  |  S(r')S(r  -r')dr' 

(B-3) 

Then  we  can  write 

y/(r)  =  G(r,r')S(r')dr' 

(B.4) 

To  find  a  solution  of  equation  (B.2),  solve  it  in  spherical  coordinates  with  the  origin  at  r' 

(V2  +  k2)G{r)  =  -S(r)  =  -S(x)S(y)S(z)  (B.5) 


Due  to  the  spherical  symmetry  of  a  point  source,  G(r )  is  also  spherical  symmetric. 

In  spherical  coordinates  equation  (B.2)  becomes  as 

-- %(rG)  +  k2G  =  -S(r)  (B.6) 

r  dr 

Everywhere,  except  at  r  =  0 ,  the  equation  is  the  same  as  the  homogeneous  equation 

^L(rG)  +  k2(rG)  =  0  (B.7) 

dr 

So  the  solution  is  shown  below 

ikr 

G(r)  =  A  —  +  B -  (B.8) 

r  r 

Since  sources  are  absent  at  infinity,  only  an  outgoing  solution  can  exist. 

ikr 

G(r)  =  A- —  (B.9) 

r 

To  find  the  constant  A  ,  substitute  equation  (B.9)  into  equation  (B.2)  and  integrate 
equation  (B.2)  over  a  small  volume. 
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1  v  ■  v 


r  r 

'  nv 


(B-10) 


Since  dV  =  4 nr1  dr  ,  the  second  integral  vanishes  and  the  first  integral  is  converted  into  a 
surface  integral  using  Gauss’  Theorem 


„  2  d  Aeikr  , 

hm4;rr" - =  -1 

r->0  dr  r 


(B  -11) 


Finally  we  get  ^4  =  V 


Therefore,  the  Greens’  function  G  is 


G(r  -  r')  = 


An\r-r'\ 


(B-12) 


Consequently,  the  final  solution  to  equation  (B.l)  is 


iK\r— r  | 

w(r)  =  ~  f  r')dr 

J  4 n \r  -  r 


(B-13) 
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APPENDIX  C.  FOURIER  TRANSFORM  IN  2-D 


When  dealing  with  both  linear  and  nonlinear  systems,  it  is  useful  to  break  down  a 
complicated  input  into  more  simple  inputs.  This  process  is  well  known  as  the  Fourier 
Transform.  The  Fourier  Transform  of  a  complex-valued  function  g  of  two  independent 
variables  x  and  y  is  defined  by 

oo 

?  {&}  =  J  J  g(x >  yy2*(V*X+Vyy)dxdy 


Similarly,  the  Inverse  Fourier  Transform  of  a  G(vx,  v  )  is  defined  as 

oo 

T-'  G(v%,vy dvxdvy  (C.2) 

-00 

But  all  functions  cannot  be  transformed  from  the  spatial  domain  into  Fourier  domain.  In 
order  to  be  transformed,  the  function  should  satisfy  “existence  conditions”  [20]. 

1 .  ‘  g  ’  must  be  absolutely  integrable  over  the  infinite  (x,  y)  plane. 

2.  ‘  g  ’  must  have  a  finite  number  of  discontinuities  and  a  finite  number  of  maxima 
and  minima  in  any  finite  rectangle. 

3.  ‘  g  ’  must  have  no  infinite  discontinuities. 

There  are  a  few  useful  properties  of  the  Fourier  Transform.  These  will  provide  the  basic 
idea  for  the  manipulation  of  the  Fourier  Transforms  and  can  save  time  solving  Fourier 
problems. 

1.  Linear  Theorem  : 

Jr(ag  +  /3k)  =  aJr(g)  +  /3Jr(k)  (C.3) 

2,  Similarity  Theorem 

If  F{g(x,y)}  =  G(vx,vy)  ,  then  F{g{ax,by)}  =  (C.4) 

\ab\  a  b 
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3. 


Shift  Theorem  : 


If  T {g(x,  v)}  =  G( vx,vy) ,  then  T [g(x-a, y-b)}  =  G{v x,v y)e'2,c(Vxa+Vyb)  (C.5) 

4.  Parseval’s  Theorem  : 


If  T  {g(x,y)}  =  G(vx,vy ) ,  then  J  J  |g(x,  y)\2dxdy  =  J  J  |  G{vx,vyf\dvxdvy  (C.6) 


5.  Convolution  Theorem  : 

If  F{g(x,y)}  =  G(vx,vy)  and  F {K(x,y)}  =  KG(vx,vy) ,  then 


(C.7) 


6.  Autocorrelation  Theorem 

If  T  {g(x,  v)}  =  G(vx,vy) ,  then 


J  J  g(a,b)g  *(a-x,b-y)dadbl  =  \G(vx,vy)\~ 


similarly,  T7  j|g(x,y)|~  j  =  J  j  \G(a,b)2\G*(a -vx,b- vy)dadb 

-oo 


(C-8) 

(C.9) 
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APPENDIX  D.  MATLAB  CODES 


%%%%%%%%%%%%  Matlab  code  for  material  reflection  %%%%%%%%%%%%% 
clear; 

%  Calculating  the  Reflection  for  plane  waves  in  various  surface 
%  Equations  used  from  Chapter  5 

%  Creating  various  index  of  refraction  for  different  materials 
%  Polystyrene  foam,  Quartz,  Glass 
nm=[1.01  1.5  1.8]; 

%  Lophira  alata(wood),  Tile,  Sand 
%  nm=[1.49  1.6  1.67]; 

%  Skin,  Body  fluid 
%  nm=[1.4  3.9]; 

%  Thickness  of  the  material 
L=0.5; 

c=3.0e8;  %  Speed  of  light 

figure(l) 

nu=linspace(lel  1 ,5e  12, 1 00); 
k=2*pi*nu./c; 

num=(  1  -nm(  1  )A2)+(nm(  1  )A2-  l)*exp(i*2*k*L*nm(  1 )); 
den=(nm(  1 )+ 1  )A2-((nm(  1 )- 1  )A2)*exp(i*2*k*L*nm(  1)); 
r=num./den; 

R=r.*conj'(r); 
hold  on 
plot(nu,R,'-'); 

gtext(sprintf(’n=%2.2f  ,nm(  1 ))); 

num=(l-nm(2)A2)+(nm(2)A2-l)*exp(i*2*k*L*nm(2)); 

den=(nm(2)+l)A2-((nm(2)-l)A2)*exp(i*2*k*L*nm(2)); 

r=num./den; 

R=r.*conj(r); 
hold  on 
plot(nu,R, 

gtext(sprintf(’n=%2 .2  f  ,nm(2 ))) ; 

num=(l-nm(3)A2)+(nm(3)A2-l)*exp(i*2*k*L*nm(3)); 

den=(nm(3)+l)A2-((nm(3)-l)A2)*exp(i*2*k*L*nm(3)); 

r=num./den; 

R=r.*conj'(r); 
hold  on 
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plot(nu,R,':'); 

gtext(sprintf('n=%2 .2f  ,nm(  3  ))) ; 
hold  off 

xlabel('frequency  (\nu)  Hz'); 
ylabel('Reflection'); 

title('Frequency  vs  Reflection  from  the  dielectric') 
gtext(sprintf(’L  =  %1.2f,L)); 

%%%%%%%%%%%  Matlab  code  for  Bom  and  Rytov  approximations  %%%%%%%% 
clear; 

%parameters 
nu=940e9; 
lambda=3 .0e8/nu; 

%k=2*  pi/lambda ; 
k=3000; 
n=1.67; 

L=0.5; 

Z=2; 

A=l; 
signal=0; 
noise=0; 

%  Defining  Coordinate  Systems  in  the  spatial  frequency  domain. 
ii=10;  %  Spatial  frequencies  for  meshgrid  input. 

fb=ii;  %  Frequency  band. 

%fn=fb; 

fn=2*fb;  %  Nyquist  rate. 

j=l/fn;  %  Separation  must  be  at  least  this  distance  to  prevent  aliasing. 

[fx,fy]=meshgrid(-ii:j:ii,-ii:j:ii);  %  Creating  coordinate  system. 

%  Defining  the  Green's  function. 

x0=0;  y0=0;  z0=-3;  %  Set  this  position  accordingly  to  the  problem. 

N=sqrt(kA2-(2*pi*fx).A2-(2*pi*fy).A2);  %  Look  in  chapter  4  for  reference. 

Nm=sqrt((n*k)A2-(2*pi*fx).A2-(2*pi*fy).A2);  %  Look  in  chapter  4  for  reference. 
num=4*N.*Nm*exp(-i*N*L);  %  Numerator  of  the  equation. 

den=(Nm+N).A2.*exp(-i*Nm*L)-(Nm-N).A2.*exp(i*Nm*L);  %  Denominator  of  the 

equation. 

%  Required  terms  for  taking  Fourier  transform 
A0=exp(-i*2*pi*fx*x0).*exp(-i*2*pi*fy*y0); 

G=(num./den).*exp(i*N*Z).*exp(i*N*z0)./(2*N).*A0;  %  Transmission  Green's  function. 
%%  Synthetic  Data  Generation 


%  Frequency  of  the  wave  propagation. 

%  In  terms  of  wave  length. 

%  Magnitude  of  the  propagation  vector. 

%  Material  index  of  refraction. 

%  Thickness  of  the  wall. 

%  Measurement  plane  at  Z=constant.  Adjust  accordingly 
%  depending  on  position  of  the  receiver. 

%  Amplitude  of  the  incident  plane  wave. 


74 


%  Transmittance  for  a  plane  wave  through  a  given  media 
t=4*n*eXp(-i*k*L)./((n+l)A2*eXp(-i*n*k*L)-(n-l)A2*eXp(i*n*k*L)); 

%  scattering  Points  in  the  spatial  frequency  domain. 

xl=-l;  yl=-l;  zl=-l; 

x2=-l;  y2=l;  z2=-3; 

x3=0;  y3=0;  z3=-5; 

x4=l;  y4=l;  z4=-1.3; 

x5=l;  y5=-l;  z5=-2.2; 

%  Incident  field  after  the  penetration  of  the  media. 

El=t*A*exp(-i*k*zl); 

E2=t*A*exp(-i*k*z2); 

E3=t*A*exp(-i*k*z3); 

E4=t*A*exp(-i*k*z4); 

E5=t*A*exp(-i*k*z5); 

%  Position  of  the  scatter  points  in  frequency  domain  for  the  given 
%  measurement  plane  at  Z=constant 

%  Since  the  scatter  points  are  in  tenns  of  dirac  delta  functions, 

%  the  field  measurement  would  become  the  summation  of  the  Green's  function 

%  and  the  incident  wave  after  performing  the  Rytov  Approximation 

Gl=(num./den).*exp(i*N*Z).*exp(i*N*zl)./(2*N); 

G2=(num./den).*exp(i*N*Z).*exp(i*N*z2)./(2*N); 

G3=(num./den).*exp(i*N*Z).*exp(i*N*z3)./(2*N); 

G4=(num./den).*exp(i*N*Z).*exp(i*N*z4)./(2*N); 

G5=(num./den).*exp(i*N*Z).*exp(i*N*z5)./(2*N); 

%  Required  terms  for  taking  Fourier  transform 

Al=exp(-i*2*pi*fx*xl).*exp(-i*2*pi*fy*yl); 

A2=exp(-i*2*pi*fx*x2).*exp(-i*2*pi*fy*y2); 

A3=exp(-i*2*pi*fx*x3).*exp(-i*2*pi*fy*y3); 

A4=exp(-i*2*pi*fx*x4).*exp(-i*2*pi*fy*y4); 

A5=exp(-i*2*pi*fx*x5).*exp(-i*2*pi*fy*y5); 

%%  Generating  Data  with  Born  approximation 

D=(E1.*G1.*A1+E2.*G2.*A2+E3.*G3.*A3+E4.*G4.*A4+E5.*G5.*A5); 

%%  Generating  Data  with  Rytov  approximation 

DD=E  1 .  *(exp(G  l.*Al)-l  )+E2.  *(exp(G2 .  *  A2)- 1  )+E3 .  *(exp(G3 .  *  A3  )- 

1  )+E4 .  *  (exp(G4 .  *  A4)- 1  )+E5 .  *(exp(G5 .  *  A5)- 1 ); 

%%  Generating  differences  by  subtracting  Born  with  Rytov 
%  Generating  Data  with  Bom 
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B=(E  1 .  *  G 1 +E2 .  *  G2+E3 .  *  G3+E4 .  *  G4+E5 .  *  G5); 

%  Generating  Data  with  Rytov 

R=E  1 .  *  (exp(G  1 )- 1  )+E2 .  *  (exp(G2)- 1  )+E3 .  *  (exp(G3 )- 1  )+E4 .  *  (exp(G4)- 
l)+E5.*(exp(G5)-l); 

Ext=(B-R).*(Al+A2+A3+A4+A5);  %  Data  with  extra  terms 

%  Setting  matrix  dimension 
[MM,NN]=size(D); 

D=D(  1  :MM,  1  :NN);  %  Change  matrix  dimension  if  warranted. 

%  Converting  bin  numbers  to  real  axis  numbers  in  spatial  domain(meters). 

%  Range  has  to  be  up  to  Nyquist  frequency. 
ax=linspace(-fn/2,fn/2,length(D)); 

%%  End  of  Data  Generation 

figure(l) 

dl=ifft2(D);  %  Data  in  the  spatial  domain. 

dl=ifftshift(dl); 

mesh(ax,ax,abs(dl));  %  Amplitude  of  field  scattered  from  5  points  using  Born 
view(-70,20) 

figure(2) 

d2=ifft2(DD);  %  Data  in  the  spatial  domain. 

d2=ifftshift(d2); 

mesh(ax,ax,abs(d2));  %  Amplitude  of  field  scattered  from  5  points  using  Rtov 
view(-70,20) 

figure(3) 

d3=ifft2(Ext);  %  Data  in  the  spatial  domain. 
d3=ifftshift(d3); 

mesh(ax,ax,abs(d3));  %  Differences  between  Born  and  Rytov  approximations 
view(-70,20) 

%%%%%%%%%%%  Matlab  code  for  regularization  mathods  %%%%%%%%%%%% 
clear; 

%parameters 

nu=940e9;  %  Frequency  of  the  wave  propagation. 

Iambda=3.0e8/nu;  %  In  terms  of  wave  length. 

%k=2*pi/lambda;  %  Magnitude  of  the  propagation  vector. 
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k=3000; 

n=1.67;  %  Material  index  of  refraction. 

L=0.5;  %  Thickness  of  the  wall. 

Z=2;  %  Measurement  plane  at  Z=constant.  Adjust  accordingly 

%  depending  on  position  of  the  receiver. 

A=1 ;  %  Amplitude  of  the  incident  plane  wave. 

signal=0; 

noise=0; 

%  Defining  Coordinate  Systems  in  the  spatial  frequency  domain. 
ii=  10;  %  Spatial  frequencies  for  meshgrid  input. 

fb=ii;  %  Frequency  band. 

fn=2*fb;  %  Nyquist  rate. 

j=l/fn;  %  Separation  must  be  at  least  this  distance  to  prevent  aliasing. 

[fx,fy]=meshgrid(-ii:j:ii,-ii:j:ii);  %  Creating  coordinate  system. 

%  Defining  the  Green's  function. 

x0=0;  y0=0;  z0=-3;  %  Set  this  position  accordingly  to  the  problem. 

N=sqrt(kA2-(2*pi*fx).A2-(2*pi*fy).A2);  %  Look  in  chapter  4  for  reference. 

Nm=sqrt((n*k)A2-(2*pi*fx).A2-(2*pi*fy).A2);  %  Look  in  chapter  4  for  reference. 

num=4*N.*Nm*exp(-i*N*L);  %  Numerator  of  the  equation. 

den=(Nm+N).A2.*exp(-i*Nm*L)-(Nm-N).A2.*exp(i*Nm*L);  %  Denominator  of  the 

equation. 

%  Required  terms  for  taking  Fourier  transform 
A0=exp(-i*2*pi*fx*x0).*exp(-i*2*pi*fy*y0); 

G=(num./den).*exp(i*N*Z).*exp(i*N*z0)./(2*N).*A0;  %  Transmission  Green's  function. 
%  Synthetic  Data  Generation 

%  Transmittance  for  a  plane  wave  through  a  given  media 
t=4*n*exp(_i*k*L)  ./((n+l)A2*exp(-i*n*k*L)-(n-l)A2*exp(i*n*k*L)); 

%  scattering  Points  in  the  spatial  frequency  domain. 

xl=-l;  yl=-l;  zl=-l; 

x2=-l;  y2=l;  z2=-3; 

x3=0;  y3=0;  z3=-5; 

x4=l;  y4=l;  z4=-1.3; 

x5=l;  y5=-l;  z5=-2.2; 

%  Incident  field  after  the  penetration  of  the  media. 

El=t*A*exp(-i*k*zl); 

E2=t*A*exp(-i*k*z2); 

E3=t*A*exp(-i*k*z3); 

E4=t*A*exp(-i*k*z4); 

E5=t*A*exp(-i*k*z5); 
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%  Position  of  the  scatter  points  in  frequency  domain  for  the  given 
%  measurement  plane  at  Z=constant 

%  Since  the  scatter  points  are  in  tenns  of  dirac  delta  functions, 

%  the  field  measurement  would  become  the  summation  of  the  Green's  function 

%  and  the  incident  wave  after  performing  the  Rytov  Approximation 

Gl=(num./den).*exp(i*N*Z).*exp(i*N*zl)./(2*N); 

G2=(num./den).*exp(i*N*Z).*exp(i*N*z2)./(2*N); 

G3=(num./den).*exp(i*N*Z).*exp(i*N*z3)./(2*N); 

G4=(num./den).*exp(i*N*Z).*exp(i*N*z4)./(2*N); 

G5=(num./den).*exp(i*N*Z).*exp(i*N*z5)./(2*N); 

%  Required  terms  for  taking  Fourier  transform 

Al=exp(-i*2*pi*fx*xl).*exp(-i*2*pi*fy*yl); 

A2=exp(-i*2*pi*fx*x2).*exp(-i*2*pi*fy*y2); 

A3=exp(-i*2*pi*fx*x3).*exp(-i*2*pi*fy*y3); 

A4=exp(-i*2*pi*fx*x4).*exp(-i*2*pi*fy*y4); 

A5=exp(-i*2*pi*fx*x5).*exp(-i*2*pi*fy*y5); 

%  Generating  Data 

D=E  1 .  *(exp(G  1 .  *  A 1 )- 1  )+E2 .  *  (exp(G2 .  *  A2)- 1  )+E3 .  *  (exp(G3 .  *  A3 )- 
1  )+E4.*(exp(G4.*A4)- 1  )+E5 .  *(exp(G5 .  *  A5)- 1 ); 

%  Setting  matrix  dimension 
[MM,NN]=size(D); 

D=D(  1  :MM,  1  :NN);  %  Change  matrix  dimension  if  warranted. 

%  Converting  bin  numbers  to  real  axis  numbers  in  spatial  domain(meters). 

%  Range  has  to  be  up  to  Nyquist  frequency. 
ax=linspace(-fn/2,fn/2,length(D)); 

n=(randn(MM,NN))+i*(randn(MM,NN));  %  noise  in  complex  values  in  the  frequency 
domain 

%  set  the  signal  to  noise  in  dB 
SNR=-20; 

S=abs(D.*conj(D)); 

NO=abs(n.*conj(n)); 
signal=sum(sum(S)); 
noise=sum(sum(NO)); 
factor=sqrt((signal/noise)*  1 0A(-SNR/ 10)); 
n=n*  factor; 

d=D+n;  %  data  with  noise  in  the  frequency  domain 

%  End  of  (noisy)  Data  Generation 

figure(2)  %  Data  corrupted  with  noise  in  spatial  domain. 
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dl=ifft2(d); 

dl=ifftshift(dl); 


%  Data  in  the  spatial  domain. 


imagesc(ax,ax,abs(dl)); 
colonnap  (1-gray); 

%grid  on 
axis([-2  2  -2  2]) 

title(sprintf('Noisy  Data  with  SNR=  %3.  If  dB',SNR)) 

%  Tikhonov  Regularization 
[u,s,v]=svd(G); 

alpha=4.6;  %  Regularization  parameter 

Reg=inv(alpha*eye(size(G'*G))+G'*G)*G'; 

figure(3) 

dd=Reg.*d; 

ddl=ifft2(dd);  %  Data  in  the  spatial  domain 

dd  1  =ifftshift(dd  1 ); 
imagesc(ax,ax,abs(dd  1 )); 
colonnap  (1-gray); 

%grid  on 
axis([-2  2-2  2]) 

title(sprintf('Tikhonov  Regularization  alpha=  %4.1f  and  SNR=%3.1f  dB', alpha, SNR)) 
%  Truncated  SVD  Method 

kappa=0.34;  %  Parameter  for  truncation  of  singular  values 

[u,s,v]=svd(G); 

%  inverse  of  the  singular  value  matrix 
sinvlu=zeros(NN,MM); 
for  j=l:rank(s) 
if  1  / s (j , j )  <  kappa 
sin  v  1  u(j  ,j  )=  1  /s  (j  ,j ) ; 
else 

sinvlu(j,j)=0; 

end 

end 


td=(v*sinvlu*u').*d; 

TD=ifft2(td);  %  Data  in  the  spatial  domain 

TD=ifftshift(TD); 


figure(4) 

imagesc(ax,ax,abs(TD)); 
colonnap  (1-gray); 


79 


axis([-2  2  -2  2]) 

title(sprintf('TSVD  w/  truncation  atkappa=  %3.2f\n  SNR=%3.1f  dB', kappa, SNR)) 
%%%%%%%%%%%%%%  Matlab  code  for  L-curve  plot  %%%%%%%%%%%%%% 


clear; 

%parameters 

nu=940e9; 

lambda=3.0e8/nu; 

%k=2*pi/lambda; 

k=3000; 

n=1.67; 

L=0.5; 

Z=2; 

A=l; 

signal=0; 

noise=0; 


%  Frequency  of  the  wave  propagation. 

%  In  terms  of  wave  length. 

%  Magnitude  of  the  propagation  vector. 

%  Material  index  of  refraction. 

%  Thickness  of  the  wall. 

%  Measurement  plane  at  Z=constant.  Adjust  accordingly 
%  depending  on  position  of  the  receiver. 

%  Amplitude  of  the  incident  plane  wave. 


%  Defining  Coordinate  Systems  in  the  spatial  frequency  domain. 
ii=10;  %  Spatial  frequencies  for  meshgrid  input. 

fb=ii;  %  Frequency  band. 

fn=2*fb;  %  Nyquist  rate. 

j=l/fn;  %  Separation  must  be  at  least  this  distance  to  prevent  aliasing. 

[fx,fy]=meshgrid(-ii:j:ii,-ii:j:ii);  %  Creating  coordinate  system. 

%  Defining  the  Green's  function. 

x0=0;  y0=0;  z0=-3;  %  Set  this  position  accordingly  to  the  problem. 

N=sqrt(kA2-(2*pi*fx).A2-(2*pi*fy).A2);  %  Look  in  chapter  4  for  reference. 

Nm=sqrt((n*k)A2-(2*pi*fx).A2-(2*pi*fy).A2);  %  Look  in  chapter  4  for  reference. 

num=4*N.*Nm*exp(-i*N*L);  %  Numerator  of  the  equation. 

den=(Nm+N).A2.*exp(-i*Nm*L)-(Nm-N).A2.*exp(i*Nm*L);  %  Denominator  of  the 

equation. 

%  Required  terms  for  taking  Fourier  transform 
A0=exp(-i*2*pi*fx*x0).*exp(-i*2*pi*fy*y0); 

%  Synthetic  Data  Generation 

%  Transmittance  for  a  plane  wave  through  a  given  media 
t=4*n*exp(-i*k*L)./((n+ 1  )A2*exp(-i*n*k*L)-(n- 1  )A2*exp(i*n*k*L)); 


%  scattering  Points  in  the  spatial  frequency  domain. 
xl=-l;  y  1  =- 1 ;  zl=-l; 
x2=-l;  y2=l;  z2=-3; 
x3=0;  y3=0;  z3=-5; 
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x4=l;  y4=l;  z4=-1.3; 
x5=l;  y5=-l;  z5=-2.2; 


%  Incident  field  after  the  penetration  of  the  media. 
El=t*A*exp(-i*k*zl); 

E2=t*A*exp(-i*k*z2); 

E3=t*A*exp(-i*k*z3); 

E4=t*A*exp(-i*k*z4); 

E5=t*A*exp(-i*k*z5); 


%  Position  of  the  scatter  points  in  frequency  domain  for  the  given 
%  measurement  plane  at  Z=constant 

%  Since  the  scatter  points  are  in  tenns  of  dirac  delta  functions, 

%  the  field  measurement  would  become  the  summation  of  the  Green's  function 

%  and  the  incident  wave  after  performing  the  Rytov  Approximation 

Gl=(num./den).*exp(i*N*Z).*exp(i*N*zl)./(2*N); 

G2=(num./den).*exp(i*N*Z).*exp(i*N*z2)./(2*N); 

G3=(num./den).*exp(i*N*Z).*exp(i*N*z3)./(2*N); 

G4=(num./den).*exp(i*N*Z).*exp(i*N*z4)./(2*N); 

G5=(num./den).*exp(i*N*Z).*exp(i*N*z5)./(2*N); 

%  Required  terms  for  taking  Fourier  transform 

Al=exp(-i*2*pi*fx*xl).*exp(-i*2*pi*fy*yl); 

A2=exp(-i*2*pi*fx*x2).*exp(-i*2*pi*fy*y2); 

A3=exp(-i*2*pi*fx*x3).*exp(-i*2*pi*fy*y3); 

A4=exp(-i*2*pi*fx*x4).*exp(-i*2*pi*fy*y4); 

A5=exp(-i*2*pi*fx*x5).*exp(-i*2*pi*fy*y5); 


G=G  1  *A  1 +G2*  A2+G3  *  A3+G4*  A4+G5  *  A5 ; 

%  Generating  Data 

D=E  1 .  *(exp(G  1 .  *  A 1 )- 1  )+E2 .  *  (exp(G2 .  *  A2)- 1  )+E3 .  *  (exp(G3  ,*A3)- 
1  )+E4 .  *  (exp(G4 .  *  A4)- 1  )+E5 .  *(exp(G5 .  *  A5)- 1 ); 

%  Setting  matrix  dimension 
[MM,NN]=size(D); 

D=D(  1  :MM,  1  :NN);  %  Change  matrix  dimension  if  warranted. 

%  Converting  bin  numbers  to  real  axis  numbers  in  spatial  domain(meters). 

%  Range  has  to  be  up  to  Nyquist  frequency. 
ax=linspace(-fn/2,fn/2,length(D)); 

no=(randn(MM,NN))+i*(randn(MM,NN));  %  noise  in  complex  values  in  the  frequency 
domain 

%  set  the  signal  to  noise  in  dB 
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SNR=-20; 

S=abs(D.*conj(D)); 

NO=abs(no.*conj(no)); 
signal=sum(sum(S)); 
noise=sum(sum(NO)); 
factor=sqrt((signal/noise)*  1 0A(-SNR/ 10)); 
no=no*  factor; 

d=D+no;  %  data  with  noise  in  the  frequency  domain 
%  End  of  (noisy)  Data  Generation 

%  Tikhonov  Regularization 
[u,s,v]=svd(G); 

for  alpha=linspace(10A-2,10A2,1000)  %  Range  of  Regularization  parameters 

Reg=inv(alpha*eye(size(G’*G))+G'*G)*G'*d;  %  Regularized  solution 
sol=norm(eye(size(G’*G))*Reg);  %  Solution  norm 

res=norm(G’*G*Reg-G’*d);  %  Residual  norm 


figure(3) 

loglog(res,sol,'-o’); 
hold  on 

title(’L-curve  for  Tikhonov  Regularization’); 

xlabel(’Residual  norm’); 

ylabel(’Solution  noun'); 

figure(4) 

plot(alpha,res); 

hold  on 

title(’alpha  VS  residual  norm’); 
xlabel(’alpha’); 
ylabel(’Residual  norm'); 

end 

title(’L-curve  for  Tikhonov  Regularization’) 
xlabel(’Residual  norm’); 
ylabel(’Solution  norm’); 
axis([10A-4  10A4  10A0  10A8]) 
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